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