source/unres/src_MIN/readrtns_min.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21 c      if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23 c      if (modecalc.eq.14) then
  24 c         call read_MDpar
  25 c         call read_REMDpar
  26 c      endif
  27 C Read MUCA control parameters if required
  28 c      if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31 c      if (modecalc.eq.8) then
  32 c       inquire (file="fort.40",exist=file_exist)
  33 c       if (.not.file_exist) call csaread
  34 c      endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39 C Print restraint information
  40 #ifdef MPI
  41       if (.not. out1file .or. me.eq.king) then
  42 #endif
  43       if (nhpb.gt.nss)
  44      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  45      & " distance constraints have been imposed"
  46       do i=nss+1,nhpb
  47         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
  48       enddo
  49 #ifdef MPI
  50       endif
  51 #endif
  52 c      print *,"Processor",myrank," leaves READRTNS"
  53       return
  54       end
  55 C-------------------------------------------------------------------------------
  56       subroutine read_control
  57 C
  58 C Read contorl data
  59 C
  60       implicit real*8 (a-h,o-z)
  61       include 'DIMENSIONS'
  62 #ifdef MP
  63       include 'mpif.h'
  64       logical OKRandom, prng_restart
  65       real*8  r1
  66 #endif
  67       include 'COMMON.IOUNITS'
  68       include 'COMMON.TIME1'
  69       include 'COMMON.SBRIDGE'
  70       include 'COMMON.CONTROL'
  71       include 'COMMON.MCM'
  72       include 'COMMON.HEADER'
  73       include 'COMMON.CHAIN'
  74       include 'COMMON.FFIELD'
  75       include 'COMMON.SETUP'
  76       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  77       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  78       character*80 ucase
  79       character*320 controlcard
  80
  81       nglob_csa=0
  82       eglob_csa=1d99
  83       nmin_csa=0
  84       read (INP,'(a)') titel
  85       call card_concat(controlcard)
  86 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  87 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  88       call reada(controlcard,'SEED',seed,0.0D0)
  89       call random_init(seed)
  90 C Set up the time limit (caution! The time must be input in minutes!)
  91       read_cart=index(controlcard,'READ_CART').gt.0
  92       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  93       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
  94       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
  95       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
  96       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
  97       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
  98       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
  99       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 100       call reada(controlcard,'DRMS',drms,0.1D0)
 101       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 102        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 103        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 104        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 105        write (iout,'(a,f10.1)')'DRMS    = ',drms
 106        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 107        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 108       endif
 109       call readi(controlcard,'NZ_START',nz_start,0)
 110       call readi(controlcard,'NZ_END',nz_end,0)
 111       call readi(controlcard,'IZ_SC',iz_sc,0)
 112       timlim=60.0D0*timlim
 113       safety = 60.0d0*safety
 114       timem=timlim
 115       modecalc=0
 116       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 117       minim=(index(controlcard,'MINIMIZE').gt.0)
 118       dccart=(index(controlcard,'CART').gt.0)
 119       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 120       overlapsc=.not.overlapsc
 121       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 122       searchsc=.not.searchsc
 123       sideadd=(index(controlcard,'SIDEADD').gt.0)
 124       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 125       outpdb=(index(controlcard,'PDBOUT').gt.0)
 126       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 127       pdbref=(index(controlcard,'PDBREF').gt.0)
 128       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 129       indpdb=index(controlcard,'PDBSTART')
 130       extconf=(index(controlcard,'EXTCONF').gt.0)
 131       call readi(controlcard,'IPRINT',iprint,0)
 132       call readi(controlcard,'MAXGEN',maxgen,10000)
 133       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 134       call readi(controlcard,"KDIAG",kdiag,0)
 135       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
 136       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 137      & write (iout,*) "RESCALE_MODE",rescale_mode
 138       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 139       if (index(controlcard,'REGULAR').gt.0.0D0) then
 140         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 141         modecalc=1
 142         refstr=.true.
 143       endif
 144       if (index(controlcard,'CHECKGRAD').gt.0) then
 145         modecalc=5
 146         if (index(controlcard,'CART').gt.0) then
 147           icheckgrad=1
 148         elseif (index(controlcard,'CARINT').gt.0) then
 149           icheckgrad=2
 150         else
 151           icheckgrad=3
 152         endif
 153       elseif (index(controlcard,'THREAD').gt.0) then
 154         modecalc=2
 155         call readi(controlcard,'THREAD',nthread,0)
 156         if (nthread.gt.0) then
 157           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 158         else
 159           if (fg_rank.eq.0)
 160      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 161           stop 'Error termination in Read_Control.'
 162         endif
 163       else if (index(controlcard,'MCMA').gt.0) then
 164         modecalc=3
 165       else if (index(controlcard,'MCEE').gt.0) then
 166         modecalc=6
 167       else if (index(controlcard,'MULTCONF').gt.0) then
 168         modecalc=4
 169       else if (index(controlcard,'MAP').gt.0) then
 170         modecalc=7
 171         call readi(controlcard,'MAP',nmap,0)
 172       else if (index(controlcard,'CSA').gt.0) then
 173         modecalc=8
 174 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 175 crc
 176 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 177 crc
 178 crc        modecalc=9
 179 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 180 cfmc        modecalc=10
 181       else if (index(controlcard,'SOFTREG').gt.0) then
 182         modecalc=11
 183       else if (index(controlcard,'CHECK_BOND').gt.0) then
 184         modecalc=-1
 185       else if (index(controlcard,'TEST').gt.0) then
 186         modecalc=-2
 187       else if (index(controlcard,'MD').gt.0) then
 188         modecalc=12
 189       else if (index(controlcard,'RE ').gt.0) then
 190         modecalc=14
 191       endif
 192
 193       lmuca=index(controlcard,'MUCA').gt.0
 194       call readi(controlcard,'MUCADYN',mucadyn,0)
 195       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 196       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 197      & then
 198        write (iout,*) 'MUCADYN=',mucadyn
 199        write (iout,*) 'MUCASMOOTH=',muca_smooth
 200       endif
 201
 202       iscode=index(controlcard,'ONE_LETTER')
 203       indphi=index(controlcard,'PHI')
 204       indback=index(controlcard,'BACK')
 205       iranconf=index(controlcard,'RAND_CONF')
 206       i2ndstr=index(controlcard,'USE_SEC_PRED')
 207       gradout=index(controlcard,'GRADOUT').gt.0
 208       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 209
 210       if(me.eq.king.or..not.out1file)
 211      & write (iout,'(2a)') diagmeth(kdiag),
 212      &  ' routine used to diagonalize matrices.'
 213       return
 214       end
 215 c------------------------------------------------------------------------------
 216       subroutine molread
 217 C
 218 C Read molecular data.
 219 C
 220       implicit real*8 (a-h,o-z)
 221       include 'DIMENSIONS'
 222 #ifdef MPI
 223       include 'mpif.h'
 224       integer error_msg
 225 #endif
 226       include 'COMMON.IOUNITS'
 227       include 'COMMON.GEO'
 228       include 'COMMON.VAR'
 229       include 'COMMON.INTERACT'
 230       include 'COMMON.LOCAL'
 231       include 'COMMON.NAMES'
 232       include 'COMMON.CHAIN'
 233       include 'COMMON.FFIELD'
 234       include 'COMMON.SBRIDGE'
 235       include 'COMMON.HEADER'
 236       include 'COMMON.CONTROL'
 237 c      include 'COMMON.DBASE'
 238 c      include 'COMMON.THREAD'
 239       include 'COMMON.CONTACTS'
 240       include 'COMMON.TORCNSTR'
 241       include 'COMMON.TIME1'
 242       include 'COMMON.BOUNDS'
 243 c      include 'COMMON.MD'
 244 c      include 'COMMON.REMD'
 245       include 'COMMON.SETUP'
 246       character*4 sequence(maxres)
 247       integer rescode
 248       double precision x(maxvar)
 249       character*256 pdbfile
 250       character*320 weightcard
 251       character*80 weightcard_t,ucase
 252       dimension itype_pdb(maxres)
 253       common /pizda/ itype_pdb
 254       logical seq_comp,fail
 255       double precision energia(0:n_ene)
 256       integer ilen
 257       external ilen
 258 C
 259 C Body
 260 C
 261 C Read weights of the subsequent energy terms.
 262
 263        call card_concat(weightcard)
 264        call reada(weightcard,'WLONG',wlong,1.0D0)
 265        call reada(weightcard,'WSC',wsc,wlong)
 266        call reada(weightcard,'WSCP',wscp,wlong)
 267        call reada(weightcard,'WELEC',welec,1.0D0)
 268        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 269        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 270        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 271        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 272        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 273        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 274        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 275        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 276        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 277        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 278        call reada(weightcard,'WBOND',wbond,1.0D0)
 279        call reada(weightcard,'WTOR',wtor,1.0D0)
 280        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 281        call reada(weightcard,'WANG',wang,1.0D0)
 282        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 283        call reada(weightcard,'SCAL14',scal14,0.4D0)
 284        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 285        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 286        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 287        call reada(weightcard,'TEMP0',temp0,300.0d0)
 288        if (index(weightcard,'SOFT').gt.0) ipot=6
 289 C 12/1/95 Added weight for the multi-body term WCORR
 290        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 291        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 292        weights(1)=wsc
 293        weights(2)=wscp
 294        weights(3)=welec
 295        weights(4)=wcorr
 296        weights(5)=wcorr5
 297        weights(6)=wcorr6
 298        weights(7)=wel_loc
 299        weights(8)=wturn3
 300        weights(9)=wturn4
 301        weights(10)=wturn6
 302        weights(11)=wang
 303        weights(12)=wscloc
 304        weights(13)=wtor
 305        weights(14)=wtor_d
 306        weights(15)=wstrain
 307        weights(16)=wvdwpp
 308        weights(17)=wbond
 309        weights(18)=scal14
 310        weights(21)=wsccor
 311
 312       if(me.eq.king.or..not.out1file)
 313      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 314      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 315      &  wturn4,wturn6
 316    10 format (/'Energy-term weights (unscaled):'//
 317      & 'WSCC=   ',f10.6,' (SC-SC)'/
 318      & 'WSCP=   ',f10.6,' (SC-p)'/
 319      & 'WELEC=  ',f10.6,' (p-p electr)'/
 320      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 321      & 'WBOND=  ',f10.6,' (stretching)'/
 322      & 'WANG=   ',f10.6,' (bending)'/
 323      & 'WSCLOC= ',f10.6,' (SC local)'/
 324      & 'WTOR=   ',f10.6,' (torsional)'/
 325      & 'WTORD=  ',f10.6,' (double torsional)'/
 326      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 327      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 328      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 329      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 330      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 331      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 332      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 333      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 334      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 335       if(me.eq.king.or..not.out1file)then
 336        if (wcorr4.gt.0.0d0) then
 337         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 338      &   'between contact pairs of peptide groups'
 339         write (iout,'(2(a,f5.3/))')
 340      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 341      &  'Range of quenching the correlation terms:',2*delt_corr
 342        else if (wcorr.gt.0.0d0) then
 343         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 344      &   'between contact pairs of peptide groups'
 345        endif
 346        write (iout,'(a,f8.3)')
 347      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 348        write (iout,'(a,f8.3)')
 349      &  'General scaling factor of SC-p interactions:',scalscp
 350       endif
 351       r0_corr=cutoff_corr-delt_corr
 352       do i=1,20
 353         aad(i,1)=scalscp*aad(i,1)
 354         aad(i,2)=scalscp*aad(i,2)
 355         bad(i,1)=scalscp*bad(i,1)
 356         bad(i,2)=scalscp*bad(i,2)
 357       enddo
 358 c      call rescale_weights(t_bath)
 359       if(me.eq.king.or..not.out1file)
 360      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 361      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 362      &  wturn4,wturn6
 363    22 format (/'Energy-term weights (scaled):'//
 364      & 'WSCC=   ',f10.6,' (SC-SC)'/
 365      & 'WSCP=   ',f10.6,' (SC-p)'/
 366      & 'WELEC=  ',f10.6,' (p-p electr)'/
 367      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 368      & 'WBOND=  ',f10.6,' (stretching)'/
 369      & 'WANG=   ',f10.6,' (bending)'/
 370      & 'WSCLOC= ',f10.6,' (SC local)'/
 371      & 'WTOR=   ',f10.6,' (torsional)'/
 372      & 'WTORD=  ',f10.6,' (double torsional)'/
 373      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 374      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 375      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 376      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 377      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 378      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 379      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 380      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 381      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 382       if(me.eq.king.or..not.out1file)
 383      & write (iout,*) "Reference temperature for weights calculation:",
 384      &  temp0
 385       call reada(weightcard,"D0CM",d0cm,3.78d0)
 386       call reada(weightcard,"AKCM",akcm,15.1d0)
 387       call reada(weightcard,"AKTH",akth,11.0d0)
 388       call reada(weightcard,"AKCT",akct,12.0d0)
 389       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 390       call reada(weightcard,"V2SS",v2ss,7.61d0)
 391       call reada(weightcard,"V3SS",v3ss,13.7d0)
 392       call reada(weightcard,"EBR",ebr,-5.50D0)
 393       if(me.eq.king.or..not.out1file) then
 394        write (iout,*) "Parameters of the SS-bond potential:"
 395        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 396      & " AKCT",akct
 397        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 398        write (iout,*) "EBR",ebr
 399        print *,'indpdb=',indpdb,' pdbref=',pdbref
 400       endif
 401       if (indpdb.gt.0 .or. pdbref) then
 402         read(inp,'(a)') pdbfile
 403         if(me.eq.king.or..not.out1file)
 404      &   write (iout,'(2a)') 'PDB data will be read from file ',
 405      &   pdbfile(:ilen(pdbfile))
 406         open(ipdbin,file=pdbfile,status='old',err=33)
 407         goto 34
 408   33    write (iout,'(a)') 'Error opening PDB file.'
 409         stop
 410   34    continue
 411 c        print *,'Begin reading pdb data'
 412         call readpdb
 413 c        print *,'Finished reading pdb data'
 414         if(me.eq.king.or..not.out1file)
 415      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 416      &   ' nstart_sup=',nstart_sup
 417         do i=1,nres
 418           itype_pdb(i)=itype(i)
 419         enddo
 420         close (ipdbin)
 421         nnt=nstart_sup
 422         nct=nstart_sup+nsup-1
 423 c        call contact(.false.,ncont_ref,icont_ref,co)
 424
 425         if (sideadd) then
 426 C Following 2 lines for diagnostics; comment out if not needed
 427          write (iout,*) "Before sideadd"
 428          call intout
 429          if(me.eq.king.or..not.out1file)
 430      &    write(iout,*)'Adding sidechains'
 431          maxsi=1000
 432          do i=2,nres-1
 433           iti=itype(i)
 434           if (iti.ne.10) then
 435             nsi=0
 436             fail=.true.
 437             do while (fail.and.nsi.le.maxsi)
 438               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 439               nsi=nsi+1
 440             enddo
 441             if(fail) write(iout,*)'Adding sidechain failed for res ',
 442      &              i,' after ',nsi,' trials'
 443           endif
 444          enddo
 445 C 9/29/12 Adam: Recalculate coordinates with new side chain positions
 446          call chainbuild
 447         endif
 448 C Following 2 lines for diagnostics; comment out if not needed
 449         write (iout,*) "After sideadd"
 450         call intout
 451       endif
 452       if (indpdb.eq.0) then
 453 C Read sequence if not taken from the pdb file.
 454         read (inp,*) nres
 455 c        print *,'nres=',nres
 456         if (iscode.gt.0) then
 457           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 458         else
 459           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 460         endif
 461 C Convert sequence to numeric code
 462         do i=1,nres
 463           itype(i)=rescode(i,sequence(i),iscode)
 464         enddo
 465 C Assign initial virtual bond lengths
 466         do i=2,nres
 467           vbld(i)=vbl
 468           vbld_inv(i)=vblinv
 469         enddo
 470         do i=2,nres-1
 471           vbld(i+nres)=dsc(itype(i))
 472           vbld_inv(i+nres)=dsc_inv(itype(i))
 473 c          write (iout,*) "i",i," itype",itype(i),
 474 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 475         enddo
 476       endif
 477 c      print *,nres
 478 c      print '(20i4)',(itype(i),i=1,nres)
 479       do i=1,nres
 480 #ifdef PROCOR
 481         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 482 #else
 483         if (itype(i).eq.21) then
 484 #endif
 485           itel(i)=0
 486 #ifdef PROCOR
 487         else if (itype(i+1).ne.20) then
 488 #else
 489         else if (itype(i).ne.20) then
 490 #endif
 491           itel(i)=1
 492         else
 493           itel(i)=2
 494         endif
 495       enddo
 496       if(me.eq.king.or..not.out1file)then
 497        write (iout,*) "ITEL"
 498        do i=1,nres-1
 499          write (iout,*) i,itype(i),itel(i)
 500        enddo
 501        print *,'Call Read_Bridge.'
 502       endif
 503       call read_bridge
 504 C 8/13/98 Set limits to generating the dihedral angles
 505       do i=1,nres
 506         phibound(1,i)=-pi
 507         phibound(2,i)=pi
 508       enddo
 509       read (inp,*) ndih_constr
 510       if (ndih_constr.gt.0) then
 511         read (inp,*) ftors
 512         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 513         if(me.eq.king.or..not.out1file)then
 514          write (iout,*)
 515      &   'There are',ndih_constr,' constraints on phi angles.'
 516          do i=1,ndih_constr
 517           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 518          enddo
 519         endif
 520         do i=1,ndih_constr
 521           phi0(i)=deg2rad*phi0(i)
 522           drange(i)=deg2rad*drange(i)
 523         enddo
 524         if(me.eq.king.or..not.out1file)
 525      &   write (iout,*) 'FTORS',ftors
 526         do i=1,ndih_constr
 527           ii = idih_constr(i)
 528           phibound(1,ii) = phi0(i)-drange(i)
 529           phibound(2,ii) = phi0(i)+drange(i)
 530         enddo
 531       endif
 532       nnt=1
 533 #ifdef MPI
 534       if (me.eq.king) then
 535 #endif
 536        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 537        do i=1,nres
 538          write (iout,'(a3,i5,2f10.1)')
 539      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 540        enddo
 541 #ifdef MP
 542       endif
 543 #endif
 544       nct=nres
 545 cd      print *,'NNT=',NNT,' NCT=',NCT
 546       if (itype(1).eq.21) nnt=2
 547       if (itype(nres).eq.21) nct=nct-1
 548       if (pdbref) then
 549         if(me.eq.king.or..not.out1file)
 550      &   write (iout,'(a,i3)') 'nsup=',nsup
 551         nstart_seq=nnt
 552         if (nsup.le.(nct-nnt+1)) then
 553           do i=0,nct-nnt+1-nsup
 554             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 555               nstart_seq=nnt+i
 556               goto 111
 557             endif
 558           enddo
 559           write (iout,'(a)')
 560      &            'Error - sequences to be superposed do not match.'
 561           stop
 562         else
 563           do i=0,nsup-(nct-nnt+1)
 564             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 565      &      then
 566               nstart_sup=nstart_sup+i
 567               nsup=nct-nnt+1
 568               goto 111
 569             endif
 570           enddo
 571           write (iout,'(a)')
 572      &            'Error - sequences to be superposed do not match.'
 573         endif
 574   111   continue
 575         if (nsup.eq.0) nsup=nct-nnt
 576         if (nstart_sup.eq.0) nstart_sup=nnt
 577         if (nstart_seq.eq.0) nstart_seq=nnt
 578         if(me.eq.king.or..not.out1file)
 579      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 580      &                 ' nstart_seq=',nstart_seq
 581       endif
 582 c--- Zscore rms -------
 583       if (nz_start.eq.0) nz_start=nnt
 584       if (nz_end.eq.0 .and. nsup.gt.0) then
 585         nz_end=nnt+nsup-1
 586       else if (nz_end.eq.0) then
 587         nz_end=nct
 588       endif
 589       if(me.eq.king.or..not.out1file)then
 590        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 591        write (iout,*) 'IZ_SC=',iz_sc
 592       endif
 593 c----------------------
 594       call init_int_table
 595       if (refstr) then
 596         if (.not.pdbref) then
 597           call read_angles(inp,*38)
 598           goto 39
 599    38     write (iout,'(a)') 'Error reading reference structure.'
 600 #ifdef MPI
 601           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 602           stop 'Error reading reference structure'
 603 #endif
 604    39     call chainbuild
 605           call setup_var
 606 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 607           nstart_sup=nnt
 608           nstart_seq=nnt
 609           nsup=nct-nnt+1
 610           do i=1,2*nres
 611             do j=1,3
 612               cref(j,i)=c(j,i)
 613             enddo
 614           enddo
 615 c          call contact(.true.,ncont_ref,icont_ref,co)
 616         endif
 617 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 618         call flush(iout)
 619         if (constr_dist.gt.0) call read_dist_constr
 620 c        write (iout,*) "After read_dist_constr nhpb",nhpb
 621         call hpb_partition
 622         if(me.eq.king.or..not.out1file)
 623      &   write (iout,*) 'Contact order:',co
 624         if (pdbref) then
 625         if(me.eq.king.or..not.out1file)
 626      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 627         do i=1,ncont_ref
 628           do j=1,2
 629             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 630           enddo
 631           if(me.eq.king.or..not.out1file)
 632      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 633      &     icont_ref(1,i),' ',
 634      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 635         enddo
 636         endif
 637       endif
 638       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 639      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 640      &    modecalc.ne.10) then
 641 C If input structure hasn't been supplied from the PDB file read or generate
 642 C initial geometry.
 643         if (iranconf.eq.0 .and. .not. extconf) then
 644           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 645      &     write (iout,'(a)') 'Initial geometry will be read in.'
 646           if (read_cart) then
 647             read(inp,'(8f10.5)',end=36,err=36)
 648      &       ((c(l,k),l=1,3),k=1,nres),
 649      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 650             call int_from_cart1(.false.)
 651             do i=1,nres-1
 652               do j=1,3
 653                 dc(j,i)=c(j,i+1)-c(j,i)
 654                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 655               enddo
 656             enddo
 657             do i=nnt,nct
 658               if (itype(i).ne.10) then
 659                 do j=1,3
 660                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 661                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 662                 enddo
 663               endif
 664             enddo
 665             return
 666           else
 667             call read_angles(inp,*36)
 668           endif
 669           goto 37
 670    36     write (iout,'(a)') 'Error reading angle file.'
 671 #ifdef MPI
 672           call mpi_finalize( MPI_COMM_WORLD,IERR )
 673 #endif
 674           stop 'Error reading angle file.'
 675    37     continue
 676         else if (extconf) then
 677          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 678      &    write (iout,'(a)') 'Extended chain initial geometry.'
 679          do i=3,nres
 680           theta(i)=90d0*deg2rad
 681          enddo
 682          do i=4,nres
 683           phi(i)=180d0*deg2rad
 684          enddo
 685          do i=2,nres-1
 686           alph(i)=110d0*deg2rad
 687          enddo
 688          do i=2,nres-1
 689           omeg(i)=-120d0*deg2rad
 690          enddo
 691         else
 692           if(me.eq.king.or..not.out1file)
 693      &     write (iout,'(a)') 'Random-generated initial geometry.'
 694
 695
 696 #ifdef MPI
 697           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 698      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 699 #endif
 700             do itrial=1,100
 701               itmp=1
 702               call gen_rand_conf(itmp,*30)
 703               goto 40
 704    30         write (iout,*) 'Failed to generate random conformation',
 705      &          ', itrial=',itrial
 706               write (*,*) 'Processor:',me,
 707      &          ' Failed to generate random conformation',
 708      &          ' itrial=',itrial
 709               call intout
 710
 711 #ifdef AIX
 712               call flush_(iout)
 713 #else
 714               call flush(iout)
 715 #endif
 716             enddo
 717             write (iout,'(a,i3,a)') 'Processor:',me,
 718      &        ' error in generating random conformation.'
 719             write (*,'(a,i3,a)') 'Processor:',me,
 720      &        ' error in generating random conformation.'
 721             call flush(iout)
 722 #ifdef MPI
 723             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 724    40       continue
 725           endif
 726 #else
 727           do itrial=1,100
 728             itmp=1
 729             call gen_rand_conf(itmp,*31)
 730             goto 40
 731    31       write (iout,*) 'Failed to generate random conformation',
 732      &        ', itrial=',itrial
 733             write (*,*) 'Failed to generate random conformation',
 734      &        ', itrial=',itrial
 735           enddo
 736           write (iout,'(a,i3,a)') 'Processor:',me,
 737      &      ' error in generating random conformation.'
 738           write (*,'(a,i3,a)') 'Processor:',me,
 739      &      ' error in generating random conformation.'
 740           stop
 741    40     continue
 742 #endif
 743         endif
 744       elseif (modecalc.eq.4) then
 745         read (inp,'(a)') intinname
 746         open (intin,file=intinname,status='old',err=333)
 747         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 748      &  write (iout,'(a)') 'intinname',intinname
 749         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 750         goto 334
 751   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 752 #ifdef MPI
 753         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 754 #endif
 755         stop 'Error opening angle file.'
 756   334   continue
 757
 758       endif
 759 C Generate distance constraints, if the PDB structure is to be regularized.
 760 c      if (nthread.gt.0) then
 761 c        call read_threadbase
 762 c      endif
 763       call setup_var
 764       if (me.eq.king .or. .not. out1file)
 765      & call intout
 766       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 767         write (iout,'(/a,i3,a)')
 768      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 769         write (iout,'(20i4)') (iss(i),i=1,ns)
 770         write (iout,'(/a/)') 'Pre-formed links are:'
 771         do i=1,nss
 772           i1=ihpb(i)-nres
 773           i2=jhpb(i)-nres
 774           it1=itype(i1)
 775           it2=itype(i2)
 776           if (me.eq.king.or..not.out1file)
 777      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 778      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 779      &    ebr,forcon(i)
 780         enddo
 781         write (iout,'(a)')
 782       endif
 783 c       if (i2ndstr.gt.0) call secstrp2dihc
 784 c      call geom_to_var(nvar,x)
 785 c      call etotal(energia(0))
 786 c      call enerprint(energia(0))
 787 c      call briefout(0,etot)
 788 c      stop
 789 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 790 cd    write (iout,'(a)') 'Variable list:'
 791 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 792 #ifdef MPI
 793       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 794      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 795      &  'Processor',myrank,': end reading molecular data.'
 796 #endif
 797       return
 798       end
 799 c--------------------------------------------------------------------------
 800       logical function seq_comp(itypea,itypeb,length)
 801       implicit none
 802       integer length,itypea(length),itypeb(length)
 803       integer i
 804       do i=1,length
 805         if (itypea(i).ne.itypeb(i)) then
 806           seq_comp=.false.
 807           return
 808         endif
 809       enddo
 810       seq_comp=.true.
 811       return
 812       end
 813 c-----------------------------------------------------------------------------
 814       subroutine read_bridge
 815 C Read information about disulfide bridges.
 816       implicit real*8 (a-h,o-z)
 817       include 'DIMENSIONS'
 818 #ifdef MPI
 819       include 'mpif.h'
 820 #endif
 821       include 'COMMON.IOUNITS'
 822       include 'COMMON.GEO'
 823       include 'COMMON.VAR'
 824       include 'COMMON.INTERACT'
 825       include 'COMMON.LOCAL'
 826       include 'COMMON.NAMES'
 827       include 'COMMON.CHAIN'
 828       include 'COMMON.FFIELD'
 829       include 'COMMON.SBRIDGE'
 830       include 'COMMON.HEADER'
 831       include 'COMMON.CONTROL'
 832 c      include 'COMMON.DBASE'
 833 c      include 'COMMON.THREAD'
 834       include 'COMMON.TIME1'
 835       include 'COMMON.SETUP'
 836 C Read bridging residues.
 837       read (inp,*) ns,(iss(i),i=1,ns)
 838       print *,'ns=',ns
 839       if(me.eq.king.or..not.out1file)
 840      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 841 C Check whether the specified bridging residues are cystines.
 842       do i=1,ns
 843         if (itype(iss(i)).ne.1) then
 844           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 845      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 846      &   ' can form a disulfide bridge?!!!'
 847           write (*,'(2a,i3,a)')
 848      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 849      &   ' can form a disulfide bridge?!!!'
 850 #ifdef MPI
 851          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 852          stop
 853 #endif
 854         endif
 855       enddo
 856 C Read preformed bridges.
 857       if (ns.gt.0) then
 858       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 859       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 860       if (nss.gt.0) then
 861         nhpb=nss
 862 C Check if the residues involved in bridges are in the specified list of
 863 C bridging residues.
 864         do i=1,nss
 865           do j=1,i-1
 866             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 867      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 868               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 869      &      ' contains residues present in other pairs.'
 870               write (*,'(a,i3,a)') 'Disulfide pair',i,
 871      &      ' contains residues present in other pairs.'
 872 #ifdef MPI
 873               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 874               stop
 875 #endif
 876             endif
 877           enddo
 878           do j=1,ns
 879             if (ihpb(i).eq.iss(j)) goto 10
 880           enddo
 881           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 882    10     continue
 883           do j=1,ns
 884             if (jhpb(i).eq.iss(j)) goto 20
 885           enddo
 886           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 887    20     continue
 888           dhpb(i)=dbr
 889           forcon(i)=fbr
 890         enddo
 891         do i=1,nss
 892           ihpb(i)=ihpb(i)+nres
 893           jhpb(i)=jhpb(i)+nres
 894         enddo
 895       endif
 896       endif
 897       return
 898       end
 899 c----------------------------------------------------------------------------
 900       subroutine read_x(kanal,*)
 901       implicit real*8 (a-h,o-z)
 902       include 'DIMENSIONS'
 903       include 'COMMON.GEO'
 904       include 'COMMON.VAR'
 905       include 'COMMON.CHAIN'
 906       include 'COMMON.IOUNITS'
 907       include 'COMMON.CONTROL'
 908       include 'COMMON.LOCAL'
 909       include 'COMMON.INTERACT'
 910 c Read coordinates from input
 911 c
 912       read(kanal,'(8f10.5)',end=10,err=10)
 913      &  ((c(l,k),l=1,3),k=1,nres),
 914      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 915       do j=1,3
 916         c(j,nres+1)=c(j,1)
 917         c(j,2*nres)=c(j,nres)
 918       enddo
 919       call int_from_cart1(.false.)
 920       do i=1,nres-1
 921         do j=1,3
 922           dc(j,i)=c(j,i+1)-c(j,i)
 923           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 924         enddo
 925       enddo
 926       do i=nnt,nct
 927         if (itype(i).ne.10) then
 928           do j=1,3
 929             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 930             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 931           enddo
 932         endif
 933       enddo
 934
 935       return
 936    10 return1
 937       end
 938 c------------------------------------------------------------------------------
 939       subroutine setup_var
 940       implicit real*8 (a-h,o-z)
 941       include 'DIMENSIONS'
 942       include 'COMMON.IOUNITS'
 943       include 'COMMON.GEO'
 944       include 'COMMON.VAR'
 945       include 'COMMON.INTERACT'
 946       include 'COMMON.LOCAL'
 947       include 'COMMON.NAMES'
 948       include 'COMMON.CHAIN'
 949       include 'COMMON.FFIELD'
 950       include 'COMMON.SBRIDGE'
 951       include 'COMMON.HEADER'
 952       include 'COMMON.CONTROL'
 953 c      include 'COMMON.DBASE'
 954 c      include 'COMMON.THREAD'
 955       include 'COMMON.TIME1'
 956 C Set up variable list.
 957       ntheta=nres-2
 958       nphi=nres-3
 959       nvar=ntheta+nphi
 960       nside=0
 961       do i=2,nres-1
 962         if (itype(i).ne.10) then
 963           nside=nside+1
 964           ialph(i,1)=nvar+nside
 965           ialph(nside,2)=i
 966         endif
 967       enddo
 968       if (indphi.gt.0) then
 969         nvar=nphi
 970       else if (indback.gt.0) then
 971         nvar=nphi+ntheta
 972       else
 973         nvar=nvar+2*nside
 974       endif
 975 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
 976       return
 977       end
 978 c----------------------------------------------------------------------------
 979       subroutine gen_dist_constr
 980 C Generate CA distance constraints.
 981       implicit real*8 (a-h,o-z)
 982       include 'DIMENSIONS'
 983       include 'COMMON.IOUNITS'
 984       include 'COMMON.GEO'
 985       include 'COMMON.VAR'
 986       include 'COMMON.INTERACT'
 987       include 'COMMON.LOCAL'
 988       include 'COMMON.NAMES'
 989       include 'COMMON.CHAIN'
 990       include 'COMMON.FFIELD'
 991       include 'COMMON.SBRIDGE'
 992       include 'COMMON.HEADER'
 993       include 'COMMON.CONTROL'
 994 c      include 'COMMON.DBASE'
 995 c      include 'COMMON.THREAD'
 996       include 'COMMON.TIME1'
 997       dimension itype_pdb(maxres)
 998       common /pizda/ itype_pdb
 999       character*2 iden
1000 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1001 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1002 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1003 cd     & ' nsup',nsup
1004       do i=nstart_sup,nstart_sup+nsup-1
1005 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1006 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1007         do j=i+2,nstart_sup+nsup-1
1008           nhpb=nhpb+1
1009           ihpb(nhpb)=i+nstart_seq-nstart_sup
1010           jhpb(nhpb)=j+nstart_seq-nstart_sup
1011           forcon(nhpb)=weidis
1012           dhpb(nhpb)=dist(i,j)
1013         enddo
1014       enddo
1015 cd      write (iout,'(a)') 'Distance constraints:'
1016 cd      do i=nss+1,nhpb
1017 cd        ii=ihpb(i)
1018 cd        jj=jhpb(i)
1019 cd        iden='CA'
1020 cd        if (ii.gt.nres) then
1021 cd          iden='SC'
1022 cd          ii=ii-nres
1023 cd          jj=jj-nres
1024 cd        endif
1025 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1026 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1027 cd     &  dhpb(i),forcon(i)
1028 cd      enddo
1029       return
1030       end
1031
1032       subroutine read_minim
1033       implicit real*8 (a-h,o-z)
1034       include 'DIMENSIONS'
1035       include 'COMMON.MINIM'
1036       include 'COMMON.IOUNITS'
1037       character*80 ucase
1038       character*320 minimcard
1039       call card_concat(minimcard)
1040       call readi(minimcard,'MAXMIN',maxmin,2000)
1041       call readi(minimcard,'MAXFUN',maxfun,5000)
1042       call readi(minimcard,'MINMIN',minmin,maxmin)
1043       call readi(minimcard,'MINFUN',minfun,maxmin)
1044       call reada(minimcard,'TOLF',tolf,1.0D-2)
1045       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1046       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1047       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1048       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1049       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1050      &         'Options in energy minimization:'
1051       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1052      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1053      & 'MinMin:',MinMin,' MinFun:',MinFun,
1054      & ' TolF:',TolF,' RTolF:',RTolF
1055       return
1056       end
1057 c----------------------------------------------------------------------------
1058       subroutine read_angles(kanal,*)
1059       implicit real*8 (a-h,o-z)
1060       include 'DIMENSIONS'
1061       include 'COMMON.GEO'
1062       include 'COMMON.VAR'
1063       include 'COMMON.CHAIN'
1064       include 'COMMON.IOUNITS'
1065       include 'COMMON.CONTROL'
1066 c Read angles from input
1067 c
1068        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1069        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1070        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1071        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1072
1073        do i=1,nres
1074 c 9/7/01 avoid 180 deg valence angle
1075         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1076 c
1077         theta(i)=deg2rad*theta(i)
1078         phi(i)=deg2rad*phi(i)
1079         alph(i)=deg2rad*alph(i)
1080         omeg(i)=deg2rad*omeg(i)
1081        enddo
1082       return
1083    10 return1
1084       end
1085 c----------------------------------------------------------------------------
1086       subroutine reada(rekord,lancuch,wartosc,default)
1087       implicit none
1088       character*(*) rekord,lancuch
1089       double precision wartosc,default
1090       integer ilen,iread
1091       external ilen
1092       iread=index(rekord,lancuch)
1093       if (iread.eq.0) then
1094         wartosc=default
1095         return
1096       endif
1097       iread=iread+ilen(lancuch)+1
1098       read (rekord(iread:),*,err=10,end=10) wartosc
1099       return
1100   10  wartosc=default
1101       return
1102       end
1103 c----------------------------------------------------------------------------
1104       subroutine readi(rekord,lancuch,wartosc,default)
1105       implicit none
1106       character*(*) rekord,lancuch
1107       integer wartosc,default
1108       integer ilen,iread
1109       external ilen
1110       iread=index(rekord,lancuch)
1111       if (iread.eq.0) then
1112         wartosc=default
1113         return
1114       endif
1115       iread=iread+ilen(lancuch)+1
1116       read (rekord(iread:),*,err=10,end=10) wartosc
1117       return
1118   10  wartosc=default
1119       return
1120       end
1121 c----------------------------------------------------------------------------
1122       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1123       implicit none
1124       integer dim,i
1125       integer tablica(dim),default
1126       character*(*) rekord,lancuch
1127       character*80 aux
1128       integer ilen,iread
1129       external ilen
1130       do i=1,dim
1131         tablica(i)=default
1132       enddo
1133       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1134       if (iread.eq.0) return
1135       iread=iread+ilen(lancuch)+1
1136       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1137    10 return
1138       end
1139 c----------------------------------------------------------------------------
1140       subroutine multreada(rekord,lancuch,tablica,dim,default)
1141       implicit none
1142       integer dim,i
1143       double precision tablica(dim),default
1144       character*(*) rekord,lancuch
1145       character*80 aux
1146       integer ilen,iread
1147       external ilen
1148       do i=1,dim
1149         tablica(i)=default
1150       enddo
1151       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1152       if (iread.eq.0) return
1153       iread=iread+ilen(lancuch)+1
1154       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1155    10 return
1156       end
1157 c----------------------------------------------------------------------------
1158       subroutine openunits
1159       implicit real*8 (a-h,o-z)
1160       include 'DIMENSIONS'
1161 #ifdef MPI
1162       include 'mpif.h'
1163       character*16 form,nodename
1164       integer nodelen
1165 #endif
1166       include 'COMMON.SETUP'
1167       include 'COMMON.IOUNITS'
1168 c      include 'COMMON.MD'
1169       include 'COMMON.CONTROL'
1170       integer lenpre,lenpot,ilen,lentmp
1171       external ilen
1172       character*3 out1file_text,ucase
1173       character*3 ll
1174       external ucase
1175 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1176       call getenv_loc("PREFIX",prefix)
1177       pref_orig = prefix
1178       call getenv_loc("POT",pot)
1179       call getenv_loc("DIRTMP",tmpdir)
1180       call getenv_loc("CURDIR",curdir)
1181       call getenv_loc("OUT1FILE",out1file_text)
1182 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1183       out1file_text=ucase(out1file_text)
1184       if (out1file_text(1:1).eq."Y") then
1185         out1file=.true.
1186       else
1187         out1file=fg_rank.gt.0
1188       endif
1189       lenpre=ilen(prefix)
1190       lenpot=ilen(pot)
1191       lentmp=ilen(tmpdir)
1192       if (lentmp.gt.0) then
1193           write (*,'(80(1h!))')
1194           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1195           write (*,'(80(1h!))')
1196           write (*,*)"All output files will be on node /tmp directory."
1197 #ifdef MPI
1198         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1199         if (me.eq.king) then
1200           write (*,*) "The master node is ",nodename
1201         else if (fg_rank.eq.0) then
1202           write (*,*) "I am the CG slave node ",nodename
1203         else
1204           write (*,*) "I am the FG slave node ",nodename
1205         endif
1206 #endif
1207         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1208         lenpre = lentmp+lenpre+1
1209       endif
1210       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1211 C Get the names and open the input files
1212 #if defined(WINIFL) || defined(WINPGI)
1213       open(1,file=pref_orig(:ilen(pref_orig))//
1214      &  '.inp',status='old',readonly,shared)
1215        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1216 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1217 C Get parameter filenames and open the parameter files.
1218       call getenv_loc('BONDPAR',bondname)
1219       open (ibond,file=bondname,status='old',readonly,shared)
1220       call getenv_loc('THETPAR',thetname)
1221       open (ithep,file=thetname,status='old',readonly,shared)
1222 #ifndef CRYST_THETA
1223       call getenv_loc('THETPARPDB',thetname_pdb)
1224       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1225 #endif
1226       call getenv_loc('ROTPAR',rotname)
1227       open (irotam,file=rotname,status='old',readonly,shared)
1228 #ifndef CRYST_SC
1229       call getenv_loc('ROTPARPDB',rotname_pdb)
1230       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1231 #endif
1232       call getenv_loc('TORPAR',torname)
1233       open (itorp,file=torname,status='old',readonly,shared)
1234       call getenv_loc('TORDPAR',tordname)
1235       open (itordp,file=tordname,status='old',readonly,shared)
1236       call getenv_loc('FOURIER',fouriername)
1237       open (ifourier,file=fouriername,status='old',readonly,shared)
1238       call getenv_loc('ELEPAR',elename)
1239       open (ielep,file=elename,status='old',readonly,shared)
1240       call getenv_loc('SIDEPAR',sidename)
1241       open (isidep,file=sidename,status='old',readonly,shared)
1242 #elif (defined CRAY) || (defined AIX)
1243       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1244      &  action='read')
1245 c      print *,"Processor",myrank," opened file 1"
1246       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1247 c      print *,"Processor",myrank," opened file 9"
1248 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1249 C Get parameter filenames and open the parameter files.
1250       call getenv_loc('BONDPAR',bondname)
1251       open (ibond,file=bondname,status='old',action='read')
1252 c      print *,"Processor",myrank," opened file IBOND"
1253       call getenv_loc('THETPAR',thetname)
1254       open (ithep,file=thetname,status='old',action='read')
1255 c      print *,"Processor",myrank," opened file ITHEP"
1256 #ifndef CRYST_THETA
1257       call getenv_loc('THETPARPDB',thetname_pdb)
1258       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1259 #endif
1260       call getenv_loc('ROTPAR',rotname)
1261       open (irotam,file=rotname,status='old',action='read')
1262 c      print *,"Processor",myrank," opened file IROTAM"
1263 #ifndef CRYST_SC
1264       call getenv_loc('ROTPARPDB',rotname_pdb)
1265       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1266 #endif
1267       call getenv_loc('TORPAR',torname)
1268       open (itorp,file=torname,status='old',action='read')
1269 c      print *,"Processor",myrank," opened file ITORP"
1270       call getenv_loc('TORDPAR',tordname)
1271       open (itordp,file=tordname,status='old',action='read')
1272 c      print *,"Processor",myrank," opened file ITORDP"
1273       call getenv_loc('SCCORPAR',sccorname)
1274       open (isccor,file=sccorname,status='old',action='read')
1275 c      print *,"Processor",myrank," opened file ISCCOR"
1276       call getenv_loc('FOURIER',fouriername)
1277       open (ifourier,file=fouriername,status='old',action='read')
1278 c      print *,"Processor",myrank," opened file IFOURIER"
1279       call getenv_loc('ELEPAR',elename)
1280       open (ielep,file=elename,status='old',action='read')
1281 c      print *,"Processor",myrank," opened file IELEP"
1282       call getenv_loc('SIDEPAR',sidename)
1283       open (isidep,file=sidename,status='old',action='read')
1284 c      print *,"Processor",myrank," opened file ISIDEP"
1285 c      print *,"Processor",myrank," opened parameter files"
1286 #elif (defined G77)
1287       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1288       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1289 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1290 C Get parameter filenames and open the parameter files.
1291       call getenv_loc('BONDPAR',bondname)
1292       open (ibond,file=bondname,status='old')
1293       call getenv_loc('THETPAR',thetname)
1294       open (ithep,file=thetname,status='old')
1295 #ifndef CRYST_THETA
1296       call getenv_loc('THETPARPDB',thetname_pdb)
1297       open (ithep_pdb,file=thetname_pdb,status='old')
1298 #endif
1299       call getenv_loc('ROTPAR',rotname)
1300       open (irotam,file=rotname,status='old')
1301 #ifndef CRYST_SC
1302       call getenv_loc('ROTPARPDB',rotname_pdb)
1303       open (irotam_pdb,file=rotname_pdb,status='old')
1304 #endif
1305       call getenv_loc('TORPAR',torname)
1306       open (itorp,file=torname,status='old')
1307       call getenv_loc('TORDPAR',tordname)
1308       open (itordp,file=tordname,status='old')
1309       call getenv_loc('SCCORPAR',sccorname)
1310       open (isccor,file=sccorname,status='old')
1311       call getenv_loc('FOURIER',fouriername)
1312       open (ifourier,file=fouriername,status='old')
1313       call getenv_loc('ELEPAR',elename)
1314       open (ielep,file=elename,status='old')
1315       call getenv_loc('SIDEPAR',sidename)
1316       open (isidep,file=sidename,status='old')
1317 #else
1318       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1319      &  readonly)
1320        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1321 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1322 C Get parameter filenames and open the parameter files.
1323       call getenv_loc('BONDPAR',bondname)
1324       open (ibond,file=bondname,status='old',readonly)
1325       call getenv_loc('THETPAR',thetname)
1326       open (ithep,file=thetname,status='old',readonly)
1327 #ifndef CRYST_THETA
1328       call getenv_loc('THETPARPDB',thetname_pdb)
1329       print *,"thetname_pdb ",thetname_pdb
1330       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1331       print *,ithep_pdb," opened"
1332 #endif
1333       call getenv_loc('ROTPAR',rotname)
1334       open (irotam,file=rotname,status='old',readonly)
1335 #ifndef CRYST_SC
1336       call getenv_loc('ROTPARPDB',rotname_pdb)
1337       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1338 #endif
1339       call getenv_loc('TORPAR',torname)
1340       open (itorp,file=torname,status='old',readonly)
1341       call getenv_loc('TORDPAR',tordname)
1342       open (itordp,file=tordname,status='old',readonly)
1343       call getenv_loc('SCCORPAR',sccorname)
1344       open (isccor,file=sccorname,status='old',readonly)
1345       call getenv_loc('FOURIER',fouriername)
1346       open (ifourier,file=fouriername,status='old',readonly)
1347       call getenv_loc('ELEPAR',elename)
1348       open (ielep,file=elename,status='old',readonly)
1349       call getenv_loc('SIDEPAR',sidename)
1350       open (isidep,file=sidename,status='old',readonly)
1351 #endif
1352 #ifndef OLDSCP
1353 C
1354 C 8/9/01 In the newest version SCp interaction constants are read from a file
1355 C Use -DOLDSCP to use hard-coded constants instead.
1356 C
1357       call getenv_loc('SCPPAR',scpname)
1358 #if defined(WINIFL) || defined(WINPGI)
1359       open (iscpp,file=scpname,status='old',readonly,shared)
1360 #elif (defined CRAY)  || (defined AIX)
1361       open (iscpp,file=scpname,status='old',action='read')
1362 #elif (defined G77)
1363       open (iscpp,file=scpname,status='old')
1364 #else
1365       open (iscpp,file=scpname,status='old',readonly)
1366 #endif
1367 #endif
1368       call getenv_loc('PATTERN',patname)
1369 #if defined(WINIFL) || defined(WINPGI)
1370       open (icbase,file=patname,status='old',readonly,shared)
1371 #elif (defined CRAY)  || (defined AIX)
1372       open (icbase,file=patname,status='old',action='read')
1373 #elif (defined G77)
1374       open (icbase,file=patname,status='old')
1375 #else
1376       open (icbase,file=patname,status='old',readonly)
1377 #endif
1378 #ifdef MPI
1379 C Open output file only for CG processes
1380 c      print *,"Processor",myrank," fg_rank",fg_rank
1381       if (fg_rank.eq.0) then
1382
1383       if (nodes.eq.1) then
1384         npos=3
1385       else
1386         npos = dlog10(dfloat(nodes-1))+1
1387       endif
1388       if (npos.lt.3) npos=3
1389       write (liczba,'(i1)') npos
1390       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1391      &  //')'
1392       write (liczba,form) me
1393       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1394      &  liczba(:ilen(liczba))
1395       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1396      &  //'.int'
1397       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1398      &  //'.pdb'
1399       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1400      &  liczba(:ilen(liczba))//'.mol2'
1401       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1402      &  liczba(:ilen(liczba))//'.stat'
1403       if (lentmp.gt.0)
1404      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1405      &      //liczba(:ilen(liczba))//'.stat')
1406       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1407      &  //'.rst'
1408 c      if(usampl) then
1409 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1410 c     & liczba(:ilen(liczba))//'.const'
1411 c      endif
1412
1413       endif
1414 #else
1415       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1416       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1417       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1418       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1419       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1420       if (lentmp.gt.0)
1421      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1422      &    '.stat')
1423       rest2name=prefix(:ilen(prefix))//'.rst'
1424 c      if(usampl) then
1425 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1426 c      endif
1427 #endif
1428 #if defined(AIX) || defined(PGI)
1429       if (me.eq.king .or. .not. out1file)
1430      &   open(iout,file=outname,status='unknown')
1431 #ifdef DEBUG
1432       if (fg_rank.gt.0) then
1433         write (liczba,'(i3.3)') myrank/nfgtasks
1434         write (ll,'(bz,i3.3)') fg_rank
1435         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1436      &   status='unknown')
1437       endif
1438 #endif
1439       if(me.eq.king) then
1440        open(igeom,file=intname,status='unknown',position='append')
1441        open(ipdb,file=pdbname,status='unknown')
1442        open(imol2,file=mol2name,status='unknown')
1443        open(istat,file=statname,status='unknown',position='append')
1444       else
1445 c1out       open(iout,file=outname,status='unknown')
1446       endif
1447 #else
1448       if (me.eq.king .or. .not.out1file)
1449      &    open(iout,file=outname,status='unknown')
1450 #ifdef DEBUG
1451       if (fg_rank.gt.0) then
1452         write (liczba,'(i3.3)') myrank/nfgtasks
1453         write (ll,'(bz,i3.3)') fg_rank
1454         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1455      &   status='unknown')
1456       endif
1457 #endif
1458       if(me.eq.king) then
1459        open(igeom,file=intname,status='unknown',access='append')
1460        open(ipdb,file=pdbname,status='unknown')
1461        open(imol2,file=mol2name,status='unknown')
1462        open(istat,file=statname,status='unknown',access='append')
1463       else
1464 c1out       open(iout,file=outname,status='unknown')
1465       endif
1466 #endif
1467       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1468       csa_seed=prefix(:lenpre)//'.CSA.seed'
1469       csa_history=prefix(:lenpre)//'.CSA.history'
1470       csa_bank=prefix(:lenpre)//'.CSA.bank'
1471       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1472       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1473       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1474 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1475       csa_int=prefix(:lenpre)//'.int'
1476       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1477       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1478       csa_in=prefix(:lenpre)//'.CSA.in'
1479 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1480 C Write file names
1481       if (me.eq.king)then
1482       write (iout,'(80(1h-))')
1483       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1484       write (iout,'(80(1h-))')
1485       write (iout,*) "Input file                      : ",
1486      &  pref_orig(:ilen(pref_orig))//'.inp'
1487       write (iout,*) "Output file                     : ",
1488      &  outname(:ilen(outname))
1489       write (iout,*)
1490       write (iout,*) "Sidechain potential file        : ",
1491      &  sidename(:ilen(sidename))
1492 #ifndef OLDSCP
1493       write (iout,*) "SCp potential file              : ",
1494      &  scpname(:ilen(scpname))
1495 #endif
1496       write (iout,*) "Electrostatic potential file    : ",
1497      &  elename(:ilen(elename))
1498       write (iout,*) "Cumulant coefficient file       : ",
1499      &  fouriername(:ilen(fouriername))
1500       write (iout,*) "Torsional parameter file        : ",
1501      &  torname(:ilen(torname))
1502       write (iout,*) "Double torsional parameter file : ",
1503      &  tordname(:ilen(tordname))
1504       write (iout,*) "SCCOR parameter file : ",
1505      &  sccorname(:ilen(sccorname))
1506       write (iout,*) "Bond & inertia constant file    : ",
1507      &  bondname(:ilen(bondname))
1508       write (iout,*) "Bending parameter file          : ",
1509      &  thetname(:ilen(thetname))
1510       write (iout,*) "Rotamer parameter file          : ",
1511      &  rotname(:ilen(rotname))
1512       write (iout,*) "Threading database              : ",
1513      &  patname(:ilen(patname))
1514       if (lentmp.ne.0)
1515      &write (iout,*)" DIRTMP                          : ",
1516      &  tmpdir(:lentmp)
1517       write (iout,'(80(1h-))')
1518       endif
1519       return
1520       end
1521 c----------------------------------------------------------------------------
1522       subroutine card_concat(card)
1523       implicit real*8 (a-h,o-z)
1524       include 'DIMENSIONS'
1525       include 'COMMON.IOUNITS'
1526       character*(*) card
1527       character*80 karta,ucase
1528       external ilen
1529       read (inp,'(a)') karta
1530       karta=ucase(karta)
1531       card=' '
1532       do while (karta(80:80).eq.'&')
1533         card=card(:ilen(card)+1)//karta(:79)
1534         read (inp,'(a)') karta
1535         karta=ucase(karta)
1536       enddo
1537       card=card(:ilen(card)+1)//karta
1538       return
1539       end
1540
1541       subroutine read_dist_constr
1542       implicit real*8 (a-h,o-z)
1543       include 'DIMENSIONS'
1544 #ifdef MPI
1545       include 'mpif.h'
1546 #endif
1547       include 'COMMON.SETUP'
1548       include 'COMMON.CONTROL'
1549       include 'COMMON.CHAIN'
1550       include 'COMMON.IOUNITS'
1551       include 'COMMON.SBRIDGE'
1552       integer ifrag_(2,100),ipair_(2,100)
1553       double precision wfrag_(100),wpair_(100)
1554       character*500 controlcard
1555 c      write (iout,*) "Calling read_dist_constr"
1556 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1557 c      call flush(iout)
1558       call card_concat(controlcard)
1559       call readi(controlcard,"NFRAG",nfrag_,0)
1560       call readi(controlcard,"NPAIR",npair_,0)
1561       call readi(controlcard,"NDIST",ndist_,0)
1562       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1563       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1564       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1565       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1566       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1567 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1568 c      write (iout,*) "IFRAG"
1569 c      do i=1,nfrag_
1570 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1571 c      enddo
1572 c      write (iout,*) "IPAIR"
1573 c      do i=1,npair_
1574 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1575 c      enddo
1576       call flush(iout)
1577       do i=1,nfrag_
1578         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1579         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1580      &    ifrag_(2,i)=nstart_sup+nsup-1
1581 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1582         call flush(iout)
1583         if (wfrag_(i).gt.0.0d0) then
1584         do j=ifrag_(1,i),ifrag_(2,i)-1
1585           do k=j+1,ifrag_(2,i)
1586             write (iout,*) "j",j," k",k
1587             ddjk=dist(j,k)
1588             if (constr_dist.eq.1) then
1589             nhpb=nhpb+1
1590             ihpb(nhpb)=j
1591             jhpb(nhpb)=k
1592               dhpb(nhpb)=ddjk
1593             forcon(nhpb)=wfrag_(i)
1594             else if (constr_dist.eq.2) then
1595               if (ddjk.le.dist_cut) then
1596                 nhpb=nhpb+1
1597                 ihpb(nhpb)=j
1598                 jhpb(nhpb)=k
1599                 dhpb(nhpb)=ddjk
1600                 forcon(nhpb)=wfrag_(i)
1601               endif
1602             else
1603               nhpb=nhpb+1
1604               ihpb(nhpb)=j
1605               jhpb(nhpb)=k
1606               dhpb(nhpb)=ddjk
1607               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1608             endif
1609 #ifdef MPI
1610             if (.not.out1file .or. me.eq.king)
1611      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1612      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1613 #else
1614             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1615      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1616 #endif
1617           enddo
1618         enddo
1619         endif
1620       enddo
1621       do i=1,npair_
1622         if (wpair_(i).gt.0.0d0) then
1623         ii = ipair_(1,i)
1624         jj = ipair_(2,i)
1625         if (ii.gt.jj) then
1626           itemp=ii
1627           ii=jj
1628           jj=itemp
1629         endif
1630         do j=ifrag_(1,ii),ifrag_(2,ii)
1631           do k=ifrag_(1,jj),ifrag_(2,jj)
1632             nhpb=nhpb+1
1633             ihpb(nhpb)=j
1634             jhpb(nhpb)=k
1635             forcon(nhpb)=wpair_(i)
1636             dhpb(nhpb)=dist(j,k)
1637 #ifdef MPI
1638             if (.not.out1file .or. me.eq.king)
1639      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1640      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1641 #else
1642             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1643      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1644 #endif
1645           enddo
1646         enddo
1647         endif
1648       enddo
1649       do i=1,ndist_
1650         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1651         if (forcon(nhpb+1).gt.0.0d0) then
1652           nhpb=nhpb+1
1653           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1654 #ifdef MPI
1655           if (.not.out1file .or. me.eq.king)
1656      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1657      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1658 #else
1659           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1660      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1661 #endif
1662         endif
1663       enddo
1664       call flush(iout)
1665       return
1666       end
1667 c-------------------------------------------------------------------------------
1668 #ifdef WINIFL
1669       subroutine flush(iu)
1670       return
1671       end
1672 #endif
1673 #ifdef AIX
1674       subroutine flush(iu)
1675       call flush_(iu)
1676       return
1677       end
1678 #endif
1679 c------------------------------------------------------------------------------
1680       subroutine copy_to_tmp(source)
1681       include "DIMENSIONS"
1682       include "COMMON.IOUNITS"
1683       character*(*) source
1684       character* 256 tmpfile
1685       integer ilen
1686       external ilen
1687       logical ex
1688       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1689       inquire(file=tmpfile,exist=ex)
1690       if (ex) then
1691         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1692      &   " to temporary directory..."
1693         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1694         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1695       endif
1696       return
1697       end
1698 c------------------------------------------------------------------------------
1699       subroutine move_from_tmp(source)
1700       include "DIMENSIONS"
1701       include "COMMON.IOUNITS"
1702       character*(*) source
1703       integer ilen
1704       external ilen
1705       write (*,*) "Moving ",source(:ilen(source)),
1706      & " from temporary directory to working directory"
1707       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1708       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1709       return
1710       end
1711 c------------------------------------------------------------------------------
1712       subroutine random_init(seed)
1713 C
1714 C Initialize random number generator
1715 C
1716       implicit real*8 (a-h,o-z)
1717       include 'DIMENSIONS'
1718 #ifdef AMD64
1719       integer*8 iseedi8
1720 #endif
1721 #ifdef MPI
1722       include 'mpif.h'
1723       logical OKRandom, prng_restart
1724       real*8  r1
1725       integer iseed_array(4)
1726 #endif
1727       include 'COMMON.IOUNITS'
1728       include 'COMMON.TIME1'
1729 c      include 'COMMON.THREAD'
1730       include 'COMMON.SBRIDGE'
1731       include 'COMMON.CONTROL'
1732       include 'COMMON.MCM'
1733 c      include 'COMMON.MAP'
1734       include 'COMMON.HEADER'
1735 c      include 'COMMON.CSA'
1736       include 'COMMON.CHAIN'
1737 c      include 'COMMON.MUCA'
1738 c      include 'COMMON.MD'
1739       include 'COMMON.FFIELD'
1740       include 'COMMON.SETUP'
1741       iseed=-dint(dabs(seed))
1742       if (iseed.eq.0) then
1743         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1744      &    'Random seed undefined. The program will stop.'
1745         write (*,'(/80(1h*)/20x,a/80(1h*))')
1746      &    'Random seed undefined. The program will stop.'
1747 #ifdef MPI
1748         call mpi_finalize(mpi_comm_world,ierr)
1749 #endif
1750         stop 'Bad random seed.'
1751       endif
1752 #ifdef MPI
1753       if (fg_rank.eq.0) then
1754       seed=seed*(me+1)+1
1755 #ifdef AMD64
1756       iseedi8=dint(seed)
1757       if(me.eq.king .or. .not. out1file)
1758      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1759       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1760       OKRandom = prng_restart(me,iseedi8)
1761 #else
1762       do i=1,4
1763        tmp=65536.0d0**(4-i)
1764        iseed_array(i) = dint(seed/tmp)
1765        seed=seed-iseed_array(i)*tmp
1766       enddo
1767       if(me.eq.king .or. .not. out1file)
1768      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1769      &                 (iseed_array(i),i=1,4)
1770       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1771      &                 (iseed_array(i),i=1,4)
1772       OKRandom = prng_restart(me,iseed_array)
1773 #endif
1774       if (OKRandom) then
1775 c        r1=ran_number(0.0D0,1.0D0)
1776         if(me.eq.king .or. .not. out1file)
1777      &   write (iout,*) 'ran_num',r1
1778         if (r1.lt.0.0d0) OKRandom=.false.
1779       endif
1780       if (.not.OKRandom) then
1781         write (iout,*) 'PRNG IS NOT WORKING!!!'
1782         print *,'PRNG IS NOT WORKING!!!'
1783         if (me.eq.0) then
1784          call flush(iout)
1785          call mpi_abort(mpi_comm_world,error_msg,ierr)
1786          stop
1787         else
1788          write (iout,*) 'too many processors for parallel prng'
1789          write (*,*) 'too many processors for parallel prng'
1790          call flush(iout)
1791          stop
1792         endif
1793       endif
1794       endif
1795 #else
1796       call vrndst(iseed)
1797 c      write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1798 #endif
1799       return
1800       end