source/unres/src_MD-M/readrtns_CSA.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       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21       if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23       if (modecalc.eq.14) then
  24          call read_MDpar
  25          call read_REMDpar
  26       endif
  27 C Read MUCA control parameters if required
  28       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       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       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.THREAD'
  70       include 'COMMON.SBRIDGE'
  71       include 'COMMON.CONTROL'
  72       include 'COMMON.MCM'
  73       include 'COMMON.MAP'
  74       include 'COMMON.HEADER'
  75       include 'COMMON.CSA'
  76       include 'COMMON.CHAIN'
  77       include 'COMMON.MUCA'
  78       include 'COMMON.MD'
  79       include 'COMMON.FFIELD'
  80       include 'COMMON.INTERACT'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call readi(controlcard,'SYM',symetr,1)
 100       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 101       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 102       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 103       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 104       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 105       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 106       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 107       call reada(controlcard,'DRMS',drms,0.1D0)
 108       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 109        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 110        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 111        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 112        write (iout,'(a,f10.1)')'DRMS    = ',drms
 113        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 114        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 115       endif
 116       call readi(controlcard,'NZ_START',nz_start,0)
 117       call readi(controlcard,'NZ_END',nz_end,0)
 118       call readi(controlcard,'IZ_SC',iz_sc,0)
 119       timlim=60.0D0*timlim
 120       safety = 60.0d0*safety
 121       timem=timlim
 122       modecalc=0
 123       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 124       minim=(index(controlcard,'MINIMIZE').gt.0)
 125       dccart=(index(controlcard,'CART').gt.0)
 126       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 127       overlapsc=.not.overlapsc
 128       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 129       searchsc=.not.searchsc
 130       sideadd=(index(controlcard,'SIDEADD').gt.0)
 131       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 132       outpdb=(index(controlcard,'PDBOUT').gt.0)
 133       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 134       pdbref=(index(controlcard,'PDBREF').gt.0)
 135       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 136       indpdb=index(controlcard,'PDBSTART')
 137       extconf=(index(controlcard,'EXTCONF').gt.0)
 138       call readi(controlcard,'IPRINT',iprint,0)
 139       call readi(controlcard,'MAXGEN',maxgen,10000)
 140       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 141       call readi(controlcard,"KDIAG",kdiag,0)
 142       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 143       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 144      & write (iout,*) "RESCALE_MODE",rescale_mode
 145       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 146       if (index(controlcard,'REGULAR').gt.0.0D0) then
 147         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 148         modecalc=1
 149         refstr=.true.
 150       endif
 151       if (index(controlcard,'CHECKGRAD').gt.0) then
 152         modecalc=5
 153         if (index(controlcard,'CART').gt.0) then
 154           icheckgrad=1
 155         elseif (index(controlcard,'CARINT').gt.0) then
 156           icheckgrad=2
 157         else
 158           icheckgrad=3
 159         endif
 160       elseif (index(controlcard,'THREAD').gt.0) then
 161         modecalc=2
 162         call readi(controlcard,'THREAD',nthread,0)
 163         if (nthread.gt.0) then
 164           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 165         else
 166           if (fg_rank.eq.0)
 167      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 168           stop 'Error termination in Read_Control.'
 169         endif
 170       else if (index(controlcard,'MCMA').gt.0) then
 171         modecalc=3
 172       else if (index(controlcard,'MCEE').gt.0) then
 173         modecalc=6
 174       else if (index(controlcard,'MULTCONF').gt.0) then
 175         modecalc=4
 176       else if (index(controlcard,'MAP').gt.0) then
 177         modecalc=7
 178         call readi(controlcard,'MAP',nmap,0)
 179       else if (index(controlcard,'CSA').gt.0) then
 180         modecalc=8
 181 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 182 crc
 183 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 184 crc
 185 crc        modecalc=9
 186 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 187 cfmc        modecalc=10
 188       else if (index(controlcard,'SOFTREG').gt.0) then
 189         modecalc=11
 190       else if (index(controlcard,'CHECK_BOND').gt.0) then
 191         modecalc=-1
 192       else if (index(controlcard,'TEST').gt.0) then
 193         modecalc=-2
 194       else if (index(controlcard,'MD').gt.0) then
 195         modecalc=12
 196       else if (index(controlcard,'RE ').gt.0) then
 197         modecalc=14
 198       endif
 199
 200       lmuca=index(controlcard,'MUCA').gt.0
 201       call readi(controlcard,'MUCADYN',mucadyn,0)
 202       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 203       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 204      & then
 205        write (iout,*) 'MUCADYN=',mucadyn
 206        write (iout,*) 'MUCASMOOTH=',muca_smooth
 207       endif
 208
 209       iscode=index(controlcard,'ONE_LETTER')
 210       indphi=index(controlcard,'PHI')
 211       indback=index(controlcard,'BACK')
 212       iranconf=index(controlcard,'RAND_CONF')
 213       i2ndstr=index(controlcard,'USE_SEC_PRED')
 214       gradout=index(controlcard,'GRADOUT').gt.0
 215       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 216       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 217       if (me.eq.king .or. .not.out1file )
 218      &  write (iout,*) "DISTCHAINMAX",distchainmax
 219
 220       if(me.eq.king.or..not.out1file)
 221      & write (iout,'(2a)') diagmeth(kdiag),
 222      &  ' routine used to diagonalize matrices.'
 223       return
 224       end
 225 c--------------------------------------------------------------------------
 226       subroutine read_REMDpar
 227 C
 228 C Read REMD settings
 229 C
 230       implicit real*8 (a-h,o-z)
 231       include 'DIMENSIONS'
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.TIME1'
 234       include 'COMMON.MD'
 235 #ifndef LANG0
 236       include 'COMMON.LANGEVIN'
 237 #else
 238       include 'COMMON.LANGEVIN.lang0'
 239 #endif
 240       include 'COMMON.INTERACT'
 241       include 'COMMON.NAMES'
 242       include 'COMMON.GEO'
 243       include 'COMMON.REMD'
 244       include 'COMMON.CONTROL'
 245       include 'COMMON.SETUP'
 246       character*80 ucase
 247       character*320 controlcard
 248       character*3200 controlcard1
 249       integer iremd_m_total
 250
 251       if(me.eq.king.or..not.out1file)
 252      & write (iout,*) "REMD setup"
 253
 254       call card_concat(controlcard)
 255       call readi(controlcard,"NREP",nrep,3)
 256       call readi(controlcard,"NSTEX",nstex,1000)
 257       call reada(controlcard,"RETMIN",retmin,10.0d0)
 258       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 259       mremdsync=(index(controlcard,'SYNC').gt.0)
 260       call readi(controlcard,"NSYN",i_sync_step,100)
 261       restart1file=(index(controlcard,'REST1FILE').gt.0)
 262       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 263       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 264       if(max_cache_traj_use.gt.max_cache_traj)
 265      &           max_cache_traj_use=max_cache_traj
 266       if(me.eq.king.or..not.out1file) then
 267 cd       if (traj1file) then
 268 crc caching is in testing - NTWX is not ignored
 269 cd        write (iout,*) "NTWX value is ignored"
 270 cd        write (iout,*) "  trajectory is stored to one file by master"
 271 cd        write (iout,*) "  before exchange at NSTEX intervals"
 272 cd       endif
 273        write (iout,*) "NREP= ",nrep
 274        write (iout,*) "NSTEX= ",nstex
 275        write (iout,*) "SYNC= ",mremdsync
 276        write (iout,*) "NSYN= ",i_sync_step
 277        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 278       endif
 279       remd_tlist=.false.
 280       if (index(controlcard,'TLIST').gt.0) then
 281          remd_tlist=.true.
 282          call card_concat(controlcard1)
 283          read(controlcard1,*) (remd_t(i),i=1,nrep)
 284          if(me.eq.king.or..not.out1file)
 285      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 286       endif
 287       remd_mlist=.false.
 288       if (index(controlcard,'MLIST').gt.0) then
 289          remd_mlist=.true.
 290          call card_concat(controlcard1)
 291          read(controlcard1,*) (remd_m(i),i=1,nrep)
 292          if(me.eq.king.or..not.out1file) then
 293           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 294           iremd_m_total=0
 295           do i=1,nrep
 296            iremd_m_total=iremd_m_total+remd_m(i)
 297           enddo
 298            write (iout,*) 'Total number of replicas ',iremd_m_total
 299           endif
 300          endif
 301       if(me.eq.king.or..not.out1file)
 302      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 303       return
 304       end
 305 c--------------------------------------------------------------------------
 306       subroutine read_MDpar
 307 C
 308 C Read MD settings
 309 C
 310       implicit real*8 (a-h,o-z)
 311       include 'DIMENSIONS'
 312       include 'COMMON.IOUNITS'
 313       include 'COMMON.TIME1'
 314       include 'COMMON.MD'
 315 #ifndef LANG0
 316       include 'COMMON.LANGEVIN'
 317 #else
 318       include 'COMMON.LANGEVIN.lang0'
 319 #endif
 320       include 'COMMON.INTERACT'
 321       include 'COMMON.NAMES'
 322       include 'COMMON.GEO'
 323       include 'COMMON.SETUP'
 324       include 'COMMON.CONTROL'
 325       include 'COMMON.SPLITELE'
 326       character*80 ucase
 327       character*320 controlcard
 328
 329       call card_concat(controlcard)
 330       call readi(controlcard,"NSTEP",n_timestep,1000000)
 331       call readi(controlcard,"NTWE",ntwe,100)
 332       call readi(controlcard,"NTWX",ntwx,1000)
 333       call reada(controlcard,"DT",d_time,1.0d-1)
 334       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 335       call reada(controlcard,"DAMAX",damax,1.0d1)
 336       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 337       call readi(controlcard,"LANG",lang,0)
 338       RESPA = index(controlcard,"RESPA") .gt. 0
 339       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 340       ntime_split0=ntime_split
 341       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 342       ntime_split0=ntime_split
 343       call reada(controlcard,"R_CUT",r_cut,2.0d0)
 344       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 345       rest = index(controlcard,"REST").gt.0
 346       tbf = index(controlcard,"TBF").gt.0
 347       usampl = index(controlcard,"USAMPL").gt.0
 348
 349       mdpdb = index(controlcard,"MDPDB").gt.0
 350       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 351       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 352       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 353       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 354       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 355       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 356       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 357       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 358       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 359       large = index(controlcard,"LARGE").gt.0
 360       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 361       rattle = index(controlcard,"RATTLE").gt.0
 362 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 363       nset=0
 364       if(usampl) then
 365         call read_fragments
 366       endif
 367
 368       if(me.eq.king.or..not.out1file) then
 369        write (iout,*)
 370        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 371        write (iout,*)
 372        write (iout,'(a)') "The units are:"
 373        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 374        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 375      &  " acceleration: angstrom/(48.9 fs)**2"
 376        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 377        write (iout,*)
 378        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 379        write (iout,'(a60,f10.5,a)')
 380      &  "Initial time step of numerical integration:",d_time,
 381      &  " natural units"
 382        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 383        if (RESPA) then
 384         write (iout,'(2a,i4,a)')
 385      &    "A-MTS algorithm used; initial time step for fast-varying",
 386      &    " short-range forces split into",ntime_split," steps."
 387         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 388      &   r_cut," lambda",rlamb
 389        endif
 390        write (iout,'(2a,f10.5)')
 391      &  "Maximum acceleration threshold to reduce the time step",
 392      &  "/increase split number:",damax
 393        write (iout,'(2a,f10.5)')
 394      &  "Maximum predicted energy drift to reduce the timestep",
 395      &  "/increase split number:",edriftmax
 396        write (iout,'(a60,f10.5)')
 397      & "Maximum velocity threshold to reduce velocities:",dvmax
 398        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 399        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 400        if (rattle) write (iout,'(a60)')
 401      &  "Rattle algorithm used to constrain the virtual bonds"
 402       endif
 403       reset_fricmat=1000
 404       if (lang.gt.0) then
 405         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 406         call reada(controlcard,"RWAT",rwat,1.4d0)
 407         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 408         surfarea=index(controlcard,"SURFAREA").gt.0
 409         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 410         if(me.eq.king.or..not.out1file)then
 411          write (iout,'(/a,$)') "Langevin dynamics calculation"
 412          if (lang.eq.1) then
 413           write (iout,'(a/)')
 414      &      " with direct integration of Langevin equations"
 415          else if (lang.eq.2) then
 416           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 417          else if (lang.eq.3) then
 418           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 419          else if (lang.eq.4) then
 420           write (iout,'(a/)') " in overdamped mode"
 421          else
 422           write (iout,'(//a,i5)')
 423      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 424           stop
 425          endif
 426          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 427          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 428          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 429          write (iout,'(a60,f10.5)')
 430      &   "Scaling factor of the friction forces:",scal_fric
 431          if (surfarea) write (iout,'(2a,i10,a)')
 432      &     "Friction coefficients will be scaled by solvent-accessible",
 433      &     " surface area every",reset_fricmat," steps."
 434         endif
 435 c Calculate friction coefficients and bounds of stochastic forces
 436         eta=6*pi*cPoise*etawat
 437         if(me.eq.king.or..not.out1file)
 438      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 439      &   ,eta
 440         gamp=scal_fric*(pstok+rwat)*eta
 441         stdfp=dsqrt(2*Rb*t_bath/d_time)
 442         do i=1,ntyp
 443           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 444           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 445         enddo
 446         if(me.eq.king.or..not.out1file)then
 447          write (iout,'(/2a/)')
 448      &   "Radii of site types and friction coefficients and std's of",
 449      &   " stochastic forces of fully exposed sites"
 450          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 451          do i=1,ntyp
 452           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 453      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 454          enddo
 455         endif
 456       else if (tbf) then
 457         if(me.eq.king.or..not.out1file)then
 458          write (iout,'(a)') "Berendsen bath calculation"
 459          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 460          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 461          if (reset_moment)
 462      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 463      &   count_reset_moment," steps"
 464          if (reset_vel)
 465      &    write (iout,'(a,i10,a)')
 466      &    "Velocities will be reset at random every",count_reset_vel,
 467      &   " steps"
 468         endif
 469       else
 470         if(me.eq.king.or..not.out1file)
 471      &   write (iout,'(a31)') "Microcanonical mode calculation"
 472       endif
 473       if(me.eq.king.or..not.out1file)then
 474        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 475        if (usampl) then
 476           write(iout,*) "MD running with constraints."
 477           write(iout,*) "Equilibration time ", eq_time, " mtus."
 478           write(iout,*) "Constraining ", nfrag," fragments."
 479           write(iout,*) "Length of each fragment, weight and q0:"
 480           do iset=1,nset
 481            write (iout,*) "Set of restraints #",iset
 482            do i=1,nfrag
 483               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 484      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 485            enddo
 486            write(iout,*) "constraints between ", npair, "fragments."
 487            write(iout,*) "constraint pairs, weights and q0:"
 488            do i=1,npair
 489             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 490      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 491            enddo
 492            write(iout,*) "angle constraints within ", nfrag_back,
 493      &      "backbone fragments."
 494            write(iout,*) "fragment, weights:"
 495            do i=1,nfrag_back
 496             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 497      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 498      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 499            enddo
 500           enddo
 501         iset=mod(kolor,nset)+1
 502        endif
 503       endif
 504       if(me.eq.king.or..not.out1file)
 505      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 506       return
 507       end
 508 c------------------------------------------------------------------------------
 509       subroutine molread
 510 C
 511 C Read molecular data.
 512 C
 513       implicit real*8 (a-h,o-z)
 514       include 'DIMENSIONS'
 515 #ifdef MPI
 516       include 'mpif.h'
 517       integer error_msg
 518 #endif
 519       include 'COMMON.IOUNITS'
 520       include 'COMMON.GEO'
 521       include 'COMMON.VAR'
 522       include 'COMMON.INTERACT'
 523       include 'COMMON.LOCAL'
 524       include 'COMMON.NAMES'
 525       include 'COMMON.CHAIN'
 526       include 'COMMON.FFIELD'
 527       include 'COMMON.SBRIDGE'
 528       include 'COMMON.HEADER'
 529       include 'COMMON.CONTROL'
 530       include 'COMMON.DBASE'
 531       include 'COMMON.THREAD'
 532       include 'COMMON.CONTACTS'
 533       include 'COMMON.TORCNSTR'
 534       include 'COMMON.TIME1'
 535       include 'COMMON.BOUNDS'
 536       include 'COMMON.MD'
 537       include 'COMMON.SETUP'
 538       character*4 sequence(maxres)
 539       integer rescode
 540       double precision x(maxvar)
 541       character*256 pdbfile
 542       character*320 weightcard
 543       character*80 weightcard_t,ucase
 544       dimension itype_pdb(maxres)
 545       common /pizda/ itype_pdb
 546       logical seq_comp,fail
 547       double precision energia(0:n_ene)
 548       integer ilen
 549       external ilen
 550 C
 551 C Body
 552 C
 553 C Read weights of the subsequent energy terms.
 554        call card_concat(weightcard)
 555        call reada(weightcard,'WLONG',wlong,1.0D0)
 556        call reada(weightcard,'WSC',wsc,wlong)
 557        call reada(weightcard,'WSCP',wscp,wlong)
 558        call reada(weightcard,'WELEC',welec,1.0D0)
 559        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 560        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 561        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 562        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 563        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 564        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 565        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 566        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 567        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 568        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 569        call reada(weightcard,'WBOND',wbond,1.0D0)
 570        call reada(weightcard,'WTOR',wtor,1.0D0)
 571        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 572        call reada(weightcard,'WANG',wang,1.0D0)
 573        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 574        call reada(weightcard,'SCAL14',scal14,0.4D0)
 575        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 576        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 577        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 578        call reada(weightcard,'TEMP0',temp0,300.0d0)
 579        if (index(weightcard,'SOFT').gt.0) ipot=6
 580 C 12/1/95 Added weight for the multi-body term WCORR
 581        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 582        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 583        weights(1)=wsc
 584        weights(2)=wscp
 585        weights(3)=welec
 586        weights(4)=wcorr
 587        weights(5)=wcorr5
 588        weights(6)=wcorr6
 589        weights(7)=wel_loc
 590        weights(8)=wturn3
 591        weights(9)=wturn4
 592        weights(10)=wturn6
 593        weights(11)=wang
 594        weights(12)=wscloc
 595        weights(13)=wtor
 596        weights(14)=wtor_d
 597        weights(15)=wstrain
 598        weights(16)=wvdwpp
 599        weights(17)=wbond
 600        weights(18)=scal14
 601        weights(21)=wsccor
 602       if(me.eq.king.or..not.out1file)
 603      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 604      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 605      &  wturn4,wturn6
 606    10 format (/'Energy-term weights (unscaled):'//
 607      & 'WSCC=   ',f10.6,' (SC-SC)'/
 608      & 'WSCP=   ',f10.6,' (SC-p)'/
 609      & 'WELEC=  ',f10.6,' (p-p electr)'/
 610      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 611      & 'WBOND=  ',f10.6,' (stretching)'/
 612      & 'WANG=   ',f10.6,' (bending)'/
 613      & 'WSCLOC= ',f10.6,' (SC local)'/
 614      & 'WTOR=   ',f10.6,' (torsional)'/
 615      & 'WTORD=  ',f10.6,' (double torsional)'/
 616      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 617      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 618      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 619      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 620      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 621      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 622      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 623      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 624      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 625       if(me.eq.king.or..not.out1file)then
 626        if (wcorr4.gt.0.0d0) then
 627         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 628      &   'between contact pairs of peptide groups'
 629         write (iout,'(2(a,f5.3/))')
 630      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 631      &  'Range of quenching the correlation terms:',2*delt_corr
 632        else if (wcorr.gt.0.0d0) then
 633         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 634      &   'between contact pairs of peptide groups'
 635        endif
 636        write (iout,'(a,f8.3)')
 637      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 638        write (iout,'(a,f8.3)')
 639      &  'General scaling factor of SC-p interactions:',scalscp
 640       endif
 641       r0_corr=cutoff_corr-delt_corr
 642       do i=1,ntyp
 643         aad(i,1)=scalscp*aad(i,1)
 644         aad(i,2)=scalscp*aad(i,2)
 645         bad(i,1)=scalscp*bad(i,1)
 646         bad(i,2)=scalscp*bad(i,2)
 647       enddo
 648       call rescale_weights(t_bath)
 649       if(me.eq.king.or..not.out1file)
 650      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 651      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 652      &  wturn4,wturn6
 653    22 format (/'Energy-term weights (scaled):'//
 654      & 'WSCC=   ',f10.6,' (SC-SC)'/
 655      & 'WSCP=   ',f10.6,' (SC-p)'/
 656      & 'WELEC=  ',f10.6,' (p-p electr)'/
 657      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 658      & 'WBOND=  ',f10.6,' (stretching)'/
 659      & 'WANG=   ',f10.6,' (bending)'/
 660      & 'WSCLOC= ',f10.6,' (SC local)'/
 661      & 'WTOR=   ',f10.6,' (torsional)'/
 662      & 'WTORD=  ',f10.6,' (double torsional)'/
 663      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 664      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 665      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 666      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 667      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 668      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 669      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 670      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 671      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 672       if(me.eq.king.or..not.out1file)
 673      & write (iout,*) "Reference temperature for weights calculation:",
 674      &  temp0
 675       call reada(weightcard,"D0CM",d0cm,3.78d0)
 676       call reada(weightcard,"AKCM",akcm,15.1d0)
 677       call reada(weightcard,"AKTH",akth,11.0d0)
 678       call reada(weightcard,"AKCT",akct,12.0d0)
 679       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 680       call reada(weightcard,"V2SS",v2ss,7.61d0)
 681       call reada(weightcard,"V3SS",v3ss,13.7d0)
 682       call reada(weightcard,"EBR",ebr,-5.50D0)
 683       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 684       do i=1,maxres
 685         dyn_ss_mask(i)=.false.
 686       enddo
 687       do i=1,maxres-1
 688         do j=i+1,maxres
 689           dyn_ssbond_ij(i,j)=1.0d300
 690         enddo
 691       enddo
 692       call reada(weightcard,"HT",Ht,0.0D0)
 693       if (dyn_ss) then
 694         ss_depth=ebr/wsc-0.25*eps(1,1)
 695         Ht=Ht/wsc-0.25*eps(1,1)
 696         akcm=akcm*wstrain/wsc
 697         akth=akth*wstrain/wsc
 698         akct=akct*wstrain/wsc
 699         v1ss=v1ss*wstrain/wsc
 700         v2ss=v2ss*wstrain/wsc
 701         v3ss=v3ss*wstrain/wsc
 702       else
 703         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 704       endif
 705
 706       if(me.eq.king.or..not.out1file) then
 707        write (iout,*) "Parameters of the SS-bond potential:"
 708        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 709      & " AKCT",akct
 710        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 711        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 712        write (iout,*)" HT",Ht
 713        print *,'indpdb=',indpdb,' pdbref=',pdbref
 714       endif
 715       if (indpdb.gt.0 .or. pdbref) then
 716         read(inp,'(a)') pdbfile
 717         if(me.eq.king.or..not.out1file)
 718      &   write (iout,'(2a)') 'PDB data will be read from file ',
 719      &   pdbfile(:ilen(pdbfile))
 720         open(ipdbin,file=pdbfile,status='old',err=33)
 721         goto 34
 722   33    write (iout,'(a)') 'Error opening PDB file.'
 723         stop
 724   34    continue
 725 c        print *,'Begin reading pdb data'
 726         call readpdb
 727 c        print *,'Finished reading pdb data'
 728         if(me.eq.king.or..not.out1file)
 729      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 730      &   ' nstart_sup=',nstart_sup
 731         do i=1,nres
 732           itype_pdb(i)=itype(i)
 733         enddo
 734         close (ipdbin)
 735         nnt=nstart_sup
 736         nct=nstart_sup+nsup-1
 737         call contact(.false.,ncont_ref,icont_ref,co)
 738
 739         if (sideadd) then
 740          if(me.eq.king.or..not.out1file)
 741      &    write(iout,*)'Adding sidechains'
 742          maxsi=1000
 743          do i=2,nres-1
 744           iti=itype(i)
 745           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 746             nsi=0
 747             fail=.true.
 748             do while (fail.and.nsi.le.maxsi)
 749               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 750               nsi=nsi+1
 751             enddo
 752             if(fail) write(iout,*)'Adding sidechain failed for res ',
 753      &              i,' after ',nsi,' trials'
 754           endif
 755          enddo
 756         endif
 757       endif
 758       if (indpdb.eq.0) then
 759 C Read sequence if not taken from the pdb file.
 760         read (inp,*) nres
 761 c        print *,'nres=',nres
 762         if (iscode.gt.0) then
 763           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 764         else
 765           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 766         endif
 767 C Convert sequence to numeric code
 768         do i=1,nres
 769           itype(i)=rescode(i,sequence(i),iscode)
 770         enddo
 771 C Assign initial virtual bond lengths
 772         do i=2,nres
 773           vbld(i)=vbl
 774           vbld_inv(i)=vblinv
 775         enddo
 776         do i=2,nres-1
 777           vbld(i+nres)=dsc(iabs(itype(i)))
 778           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 779 c          write (iout,*) "i",i," itype",itype(i),
 780 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 781         enddo
 782       endif
 783 c      print *,nres
 784 c      print '(20i4)',(itype(i),i=1,nres)
 785       do i=1,nres
 786 #ifdef PROCOR
 787         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 788 #else
 789         if (itype(i).eq.ntyp1) then
 790 #endif
 791           itel(i)=0
 792 #ifdef PROCOR
 793         else if (iabs(itype(i+1)).ne.20) then
 794 #else
 795         else if (iabs(itype(i)).ne.20) then
 796 #endif
 797           itel(i)=1
 798         else
 799           itel(i)=2
 800         endif
 801       enddo
 802       if(me.eq.king.or..not.out1file)then
 803        write (iout,*) "ITEL"
 804        do i=1,nres-1
 805          write (iout,*) i,itype(i),itel(i)
 806        enddo
 807        print *,'Call Read_Bridge.'
 808       endif
 809       call read_bridge
 810 C 8/13/98 Set limits to generating the dihedral angles
 811       do i=1,nres
 812         phibound(1,i)=-pi
 813         phibound(2,i)=pi
 814       enddo
 815       read (inp,*) ndih_constr
 816       if (ndih_constr.gt.0) then
 817         read (inp,*) ftors
 818         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 819         if(me.eq.king.or..not.out1file)then
 820          write (iout,*)
 821      &   'There are',ndih_constr,' constraints on phi angles.'
 822          do i=1,ndih_constr
 823           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 824          enddo
 825         endif
 826         do i=1,ndih_constr
 827           phi0(i)=deg2rad*phi0(i)
 828           drange(i)=deg2rad*drange(i)
 829         enddo
 830         if(me.eq.king.or..not.out1file)
 831      &   write (iout,*) 'FTORS',ftors
 832         do i=1,ndih_constr
 833           ii = idih_constr(i)
 834           phibound(1,ii) = phi0(i)-drange(i)
 835           phibound(2,ii) = phi0(i)+drange(i)
 836         enddo
 837       endif
 838       nnt=1
 839 #ifdef MPI
 840       if (me.eq.king) then
 841 #endif
 842        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 843        do i=1,nres
 844          write (iout,'(a3,i5,2f10.1)')
 845      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 846        enddo
 847 #ifdef MP
 848       endif
 849 #endif
 850       nct=nres
 851 cd      print *,'NNT=',NNT,' NCT=',NCT
 852       if (itype(1).eq.ntyp1) nnt=2
 853       if (itype(nres).eq.ntyp1) nct=nct-1
 854       if (pdbref) then
 855         if(me.eq.king.or..not.out1file)
 856      &   write (iout,'(a,i3)') 'nsup=',nsup
 857         nstart_seq=nnt
 858         if (nsup.le.(nct-nnt+1)) then
 859           do i=0,nct-nnt+1-nsup
 860             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 861               nstart_seq=nnt+i
 862               goto 111
 863             endif
 864           enddo
 865           write (iout,'(a)')
 866      &            'Error - sequences to be superposed do not match.'
 867           stop
 868         else
 869           do i=0,nsup-(nct-nnt+1)
 870             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 871      &      then
 872               nstart_sup=nstart_sup+i
 873               nsup=nct-nnt+1
 874               goto 111
 875             endif
 876           enddo
 877           write (iout,'(a)')
 878      &            'Error - sequences to be superposed do not match.'
 879         endif
 880   111   continue
 881         if (nsup.eq.0) nsup=nct-nnt
 882         if (nstart_sup.eq.0) nstart_sup=nnt
 883         if (nstart_seq.eq.0) nstart_seq=nnt
 884         if(me.eq.king.or..not.out1file)
 885      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 886      &                 ' nstart_seq=',nstart_seq
 887       endif
 888 c--- Zscore rms -------
 889       if (nz_start.eq.0) nz_start=nnt
 890       if (nz_end.eq.0 .and. nsup.gt.0) then
 891         nz_end=nnt+nsup-1
 892       else if (nz_end.eq.0) then
 893         nz_end=nct
 894       endif
 895       if(me.eq.king.or..not.out1file)then
 896        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 897        write (iout,*) 'IZ_SC=',iz_sc
 898       endif
 899 c----------------------
 900       call init_int_table
 901       if (refstr) then
 902         if (.not.pdbref) then
 903           call read_angles(inp,*38)
 904           goto 39
 905    38     write (iout,'(a)') 'Error reading reference structure.'
 906 #ifdef MPI
 907           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 908           stop 'Error reading reference structure'
 909 #endif
 910    39     call chainbuild
 911           call setup_var
 912 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 913           nstart_sup=nnt
 914           nstart_seq=nnt
 915           nsup=nct-nnt+1
 916           kkk=1
 917           do i=1,2*nres
 918             do j=1,3
 919               cref(j,i,kkk)=c(j,i)
 920             enddo
 921           enddo
 922           call contact(.true.,ncont_ref,icont_ref,co)
 923         endif
 924 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 925         call flush(iout)
 926         if (constr_dist.gt.0) call read_dist_constr
 927         write (iout,*) "After read_dist_constr nhpb",nhpb
 928         call hpb_partition
 929         if(me.eq.king.or..not.out1file)
 930      &   write (iout,*) 'Contact order:',co
 931         if (pdbref) then
 932         if(me.eq.king.or..not.out1file)
 933      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 934         do i=1,ncont_ref
 935           do j=1,2
 936             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 937           enddo
 938           if(me.eq.king.or..not.out1file)
 939      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 940      &     icont_ref(1,i),' ',
 941      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 942         enddo
 943         endif
 944       endif
 945       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 946      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 947      &    modecalc.ne.10) then
 948 C If input structure hasn't been supplied from the PDB file read or generate
 949 C initial geometry.
 950         if (iranconf.eq.0 .and. .not. extconf) then
 951           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 952      &     write (iout,'(a)') 'Initial geometry will be read in.'
 953           if (read_cart) then
 954             read(inp,'(8f10.5)',end=36,err=36)
 955      &       ((c(l,k),l=1,3),k=1,nres),
 956      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 957             write (iout,*) "Exit READ_CART"
 958             write (iout,'(8f10.5)')
 959      &       ((c(l,k),l=1,3),k=1,nres),
 960      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 961             call int_from_cart1(.true.)
 962             write (iout,*) "Finish INT_TO_CART"
 963             do i=1,nres-1
 964               do j=1,3
 965                 dc(j,i)=c(j,i+1)-c(j,i)
 966                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 967               enddo
 968             enddo
 969             do i=nnt,nct
 970               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 971                 do j=1,3
 972                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 973                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 974                 enddo
 975               endif
 976             enddo
 977             return
 978           else
 979             call read_angles(inp,*36)
 980           endif
 981           goto 37
 982    36     write (iout,'(a)') 'Error reading angle file.'
 983 #ifdef MPI
 984           call mpi_finalize( MPI_COMM_WORLD,IERR )
 985 #endif
 986           stop 'Error reading angle file.'
 987    37     continue
 988         else if (extconf) then
 989          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 990      &    write (iout,'(a)') 'Extended chain initial geometry.'
 991          do i=3,nres
 992           theta(i)=90d0*deg2rad
 993          enddo
 994          do i=4,nres
 995           phi(i)=180d0*deg2rad
 996          enddo
 997          do i=2,nres-1
 998           alph(i)=110d0*deg2rad
 999          enddo
1000          do i=2,nres-1
1001           omeg(i)=-120d0*deg2rad
1002           if (itype(i).le.0) omeg(i)=-omeg(i)
1003          enddo
1004         else
1005           if(me.eq.king.or..not.out1file)
1006      &     write (iout,'(a)') 'Random-generated initial geometry.'
1007
1008
1009 #ifdef MPI
1010           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1011      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1012 #endif
1013             do itrial=1,100
1014               itmp=1
1015               call gen_rand_conf(itmp,*30)
1016               goto 40
1017    30         write (iout,*) 'Failed to generate random conformation',
1018      &          ', itrial=',itrial
1019               write (*,*) 'Processor:',me,
1020      &          ' Failed to generate random conformation',
1021      &          ' itrial=',itrial
1022               call intout
1023
1024 #ifdef AIX
1025               call flush_(iout)
1026 #else
1027               call flush(iout)
1028 #endif
1029             enddo
1030             write (iout,'(a,i3,a)') 'Processor:',me,
1031      &        ' error in generating random conformation.'
1032             write (*,'(a,i3,a)') 'Processor:',me,
1033      &        ' error in generating random conformation.'
1034             call flush(iout)
1035 #ifdef MPI
1036             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1037    40       continue
1038           endif
1039 #else
1040           write (*,'(a)')
1041      &      ' error in generating random conformation.'
1042           stop
1043    40     continue
1044 #endif
1045         endif
1046       elseif (modecalc.eq.4) then
1047         read (inp,'(a)') intinname
1048         open (intin,file=intinname,status='old',err=333)
1049         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1050      &  write (iout,'(a)') 'intinname',intinname
1051         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1052         goto 334
1053   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1054 #ifdef MPI
1055         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1056 #endif
1057         stop 'Error opening angle file.'
1058   334   continue
1059
1060       endif
1061 C Generate distance constraints, if the PDB structure is to be regularized.
1062       if (nthread.gt.0) then
1063         call read_threadbase
1064       endif
1065       call setup_var
1066       if (me.eq.king .or. .not. out1file)
1067      & call intout
1068       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1069         write (iout,'(/a,i3,a)')
1070      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1071         write (iout,'(20i4)') (iss(i),i=1,ns)
1072        if (dyn_ss) then
1073           write(iout,*)"Running with dynamic disulfide-bond formation"
1074        else
1075         write (iout,'(/a/)') 'Pre-formed links are:'
1076         do i=1,nss
1077           i1=ihpb(i)-nres
1078           i2=jhpb(i)-nres
1079           it1=itype(i1)
1080           it2=itype(i2)
1081           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1082      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1083      &    ebr,forcon(i)
1084         enddo
1085         write (iout,'(a)')
1086        endif
1087       endif
1088       if (ns.gt.0.and.dyn_ss) then
1089           do i=nss+1,nhpb
1090             ihpb(i-nss)=ihpb(i)
1091             jhpb(i-nss)=jhpb(i)
1092             forcon(i-nss)=forcon(i)
1093             dhpb(i-nss)=dhpb(i)
1094           enddo
1095           nhpb=nhpb-nss
1096           nss=0
1097           call hpb_partition
1098           do i=1,ns
1099             dyn_ss_mask(iss(i))=.true.
1100           enddo
1101       endif
1102       if (i2ndstr.gt.0) call secstrp2dihc
1103 c      call geom_to_var(nvar,x)
1104 c      call etotal(energia(0))
1105 c      call enerprint(energia(0))
1106 c      call briefout(0,etot)
1107 c      stop
1108 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1109 cd    write (iout,'(a)') 'Variable list:'
1110 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1111 #ifdef MPI
1112       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1113      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1114      &  'Processor',myrank,': end reading molecular data.'
1115 #endif
1116       return
1117       end
1118 c--------------------------------------------------------------------------
1119       logical function seq_comp(itypea,itypeb,length)
1120       implicit none
1121       integer length,itypea(length),itypeb(length)
1122       integer i
1123       do i=1,length
1124         if (itypea(i).ne.itypeb(i)) then
1125           seq_comp=.false.
1126           return
1127         endif
1128       enddo
1129       seq_comp=.true.
1130       return
1131       end
1132 c-----------------------------------------------------------------------------
1133       subroutine read_bridge
1134 C Read information about disulfide bridges.
1135       implicit real*8 (a-h,o-z)
1136       include 'DIMENSIONS'
1137 #ifdef MPI
1138       include 'mpif.h'
1139 #endif
1140       include 'COMMON.IOUNITS'
1141       include 'COMMON.GEO'
1142       include 'COMMON.VAR'
1143       include 'COMMON.INTERACT'
1144       include 'COMMON.LOCAL'
1145       include 'COMMON.NAMES'
1146       include 'COMMON.CHAIN'
1147       include 'COMMON.FFIELD'
1148       include 'COMMON.SBRIDGE'
1149       include 'COMMON.HEADER'
1150       include 'COMMON.CONTROL'
1151       include 'COMMON.DBASE'
1152       include 'COMMON.THREAD'
1153       include 'COMMON.TIME1'
1154       include 'COMMON.SETUP'
1155 C Read bridging residues.
1156       read (inp,*) ns,(iss(i),i=1,ns)
1157       print *,'ns=',ns
1158       if(me.eq.king.or..not.out1file)
1159      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1160 C Check whether the specified bridging residues are cystines.
1161       do i=1,ns
1162         if (itype(iss(i)).ne.1) then
1163           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1164      &   'Do you REALLY think that the residue ',
1165      &    restyp(itype(iss(i))),i,
1166      &   ' can form a disulfide bridge?!!!'
1167           write (*,'(2a,i3,a)')
1168      &   'Do you REALLY think that the residue ',
1169      &    restyp(itype(iss(i))),i,
1170      &   ' can form a disulfide bridge?!!!'
1171 #ifdef MPI
1172          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1173          stop
1174 #endif
1175         endif
1176       enddo
1177 C Read preformed bridges.
1178       if (ns.gt.0) then
1179       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1180       if(fg_rank.eq.0)
1181      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1182       if (nss.gt.0) then
1183         nhpb=nss
1184 C Check if the residues involved in bridges are in the specified list of
1185 C bridging residues.
1186         do i=1,nss
1187           do j=1,i-1
1188             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1189      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1190               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1191      &      ' contains residues present in other pairs.'
1192               write (*,'(a,i3,a)') 'Disulfide pair',i,
1193      &      ' contains residues present in other pairs.'
1194 #ifdef MPI
1195               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1196               stop
1197 #endif
1198             endif
1199           enddo
1200           do j=1,ns
1201             if (ihpb(i).eq.iss(j)) goto 10
1202           enddo
1203           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1204    10     continue
1205           do j=1,ns
1206             if (jhpb(i).eq.iss(j)) goto 20
1207           enddo
1208           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1209    20     continue
1210           dhpb(i)=dbr
1211           forcon(i)=fbr
1212         enddo
1213         do i=1,nss
1214           ihpb(i)=ihpb(i)+nres
1215           jhpb(i)=jhpb(i)+nres
1216         enddo
1217       endif
1218       endif
1219       return
1220       end
1221 c----------------------------------------------------------------------------
1222       subroutine read_x(kanal,*)
1223       implicit real*8 (a-h,o-z)
1224       include 'DIMENSIONS'
1225       include 'COMMON.GEO'
1226       include 'COMMON.VAR'
1227       include 'COMMON.CHAIN'
1228       include 'COMMON.IOUNITS'
1229       include 'COMMON.CONTROL'
1230       include 'COMMON.LOCAL'
1231       include 'COMMON.INTERACT'
1232 c Read coordinates from input
1233 c
1234       read(kanal,'(8f10.5)',end=10,err=10)
1235      &  ((c(l,k),l=1,3),k=1,nres),
1236      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1237       do j=1,3
1238         c(j,nres+1)=c(j,1)
1239         c(j,2*nres)=c(j,nres)
1240       enddo
1241       call int_from_cart1(.false.)
1242       do i=1,nres-1
1243         do j=1,3
1244           dc(j,i)=c(j,i+1)-c(j,i)
1245           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1246         enddo
1247       enddo
1248       do i=nnt,nct
1249         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1250           do j=1,3
1251             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1252             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1253           enddo
1254         endif
1255       enddo
1256
1257       return
1258    10 return1
1259       end
1260 c----------------------------------------------------------------------------
1261       subroutine read_threadbase
1262       implicit real*8 (a-h,o-z)
1263       include 'DIMENSIONS'
1264       include 'COMMON.IOUNITS'
1265       include 'COMMON.GEO'
1266       include 'COMMON.VAR'
1267       include 'COMMON.INTERACT'
1268       include 'COMMON.LOCAL'
1269       include 'COMMON.NAMES'
1270       include 'COMMON.CHAIN'
1271       include 'COMMON.FFIELD'
1272       include 'COMMON.SBRIDGE'
1273       include 'COMMON.HEADER'
1274       include 'COMMON.CONTROL'
1275       include 'COMMON.DBASE'
1276       include 'COMMON.THREAD'
1277       include 'COMMON.TIME1'
1278 C Read pattern database for threading.
1279       read (icbase,*) nseq
1280       do i=1,nseq
1281         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1282      &   nres_base(2,i),nres_base(3,i)
1283         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1284      &   nres_base(1,i))
1285 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1286 c    &   nres_base(2,i),nres_base(3,i)
1287 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1288 c    &   nres_base(1,i))
1289       enddo
1290       close (icbase)
1291       if (weidis.eq.0.0D0) weidis=0.1D0
1292       do i=nnt,nct
1293         do j=i+2,nct
1294           nhpb=nhpb+1
1295           ihpb(nhpb)=i
1296           jhpb(nhpb)=j
1297           forcon(nhpb)=weidis
1298         enddo
1299       enddo
1300       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1301       write (iout,'(a,i5)') 'nexcl: ',nexcl
1302       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1303       return
1304       end
1305 c------------------------------------------------------------------------------
1306       subroutine setup_var
1307       implicit real*8 (a-h,o-z)
1308       include 'DIMENSIONS'
1309       include 'COMMON.IOUNITS'
1310       include 'COMMON.GEO'
1311       include 'COMMON.VAR'
1312       include 'COMMON.INTERACT'
1313       include 'COMMON.LOCAL'
1314       include 'COMMON.NAMES'
1315       include 'COMMON.CHAIN'
1316       include 'COMMON.FFIELD'
1317       include 'COMMON.SBRIDGE'
1318       include 'COMMON.HEADER'
1319       include 'COMMON.CONTROL'
1320       include 'COMMON.DBASE'
1321       include 'COMMON.THREAD'
1322       include 'COMMON.TIME1'
1323 C Set up variable list.
1324       ntheta=nres-2
1325       nphi=nres-3
1326       nvar=ntheta+nphi
1327       nside=0
1328       do i=2,nres-1
1329         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1330           nside=nside+1
1331           ialph(i,1)=nvar+nside
1332           ialph(nside,2)=i
1333         endif
1334       enddo
1335       if (indphi.gt.0) then
1336         nvar=nphi
1337       else if (indback.gt.0) then
1338         nvar=nphi+ntheta
1339       else
1340         nvar=nvar+2*nside
1341       endif
1342 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1343       return
1344       end
1345 c----------------------------------------------------------------------------
1346       subroutine gen_dist_constr
1347 C Generate CA distance constraints.
1348       implicit real*8 (a-h,o-z)
1349       include 'DIMENSIONS'
1350       include 'COMMON.IOUNITS'
1351       include 'COMMON.GEO'
1352       include 'COMMON.VAR'
1353       include 'COMMON.INTERACT'
1354       include 'COMMON.LOCAL'
1355       include 'COMMON.NAMES'
1356       include 'COMMON.CHAIN'
1357       include 'COMMON.FFIELD'
1358       include 'COMMON.SBRIDGE'
1359       include 'COMMON.HEADER'
1360       include 'COMMON.CONTROL'
1361       include 'COMMON.DBASE'
1362       include 'COMMON.THREAD'
1363       include 'COMMON.TIME1'
1364       dimension itype_pdb(maxres)
1365       common /pizda/ itype_pdb
1366       character*2 iden
1367 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1368 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1369 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1370 cd     & ' nsup',nsup
1371       do i=nstart_sup,nstart_sup+nsup-1
1372 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1373 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1374         do j=i+2,nstart_sup+nsup-1
1375           nhpb=nhpb+1
1376           ihpb(nhpb)=i+nstart_seq-nstart_sup
1377           jhpb(nhpb)=j+nstart_seq-nstart_sup
1378           forcon(nhpb)=weidis
1379           dhpb(nhpb)=dist(i,j)
1380         enddo
1381       enddo
1382 cd      write (iout,'(a)') 'Distance constraints:'
1383 cd      do i=nss+1,nhpb
1384 cd        ii=ihpb(i)
1385 cd        jj=jhpb(i)
1386 cd        iden='CA'
1387 cd        if (ii.gt.nres) then
1388 cd          iden='SC'
1389 cd          ii=ii-nres
1390 cd          jj=jj-nres
1391 cd        endif
1392 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1393 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1394 cd     &  dhpb(i),forcon(i)
1395 cd      enddo
1396       return
1397       end
1398 c----------------------------------------------------------------------------
1399       subroutine map_read
1400       implicit real*8 (a-h,o-z)
1401       include 'DIMENSIONS'
1402       include 'COMMON.MAP'
1403       include 'COMMON.IOUNITS'
1404       character*3 angid(4) /'THE','PHI','ALP','OME'/
1405       character*80 mapcard,ucase
1406       do imap=1,nmap
1407         read (inp,'(a)') mapcard
1408         mapcard=ucase(mapcard)
1409         if (index(mapcard,'PHI').gt.0) then
1410           kang(imap)=1
1411         else if (index(mapcard,'THE').gt.0) then
1412           kang(imap)=2
1413         else if (index(mapcard,'ALP').gt.0) then
1414           kang(imap)=3
1415         else if (index(mapcard,'OME').gt.0) then
1416           kang(imap)=4
1417         else
1418           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1419           stop 'Error - illegal variable spec in MAP card.'
1420         endif
1421         call readi (mapcard,'RES1',res1(imap),0)
1422         call readi (mapcard,'RES2',res2(imap),0)
1423         if (res1(imap).eq.0) then
1424           res1(imap)=res2(imap)
1425         else if (res2(imap).eq.0) then
1426           res2(imap)=res1(imap)
1427         endif
1428         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1429           write (iout,'(a)')
1430      &    'Error - illegal definition of variable group in MAP.'
1431           stop 'Error - illegal definition of variable group in MAP.'
1432         endif
1433         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1434         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1435         call readi(mapcard,'NSTEP',nstep(imap),0)
1436         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1437           write (iout,'(a)')
1438      &     'Illegal boundary and/or step size specification in MAP.'
1439           stop 'Illegal boundary and/or step size specification in MAP.'
1440         endif
1441       enddo ! imap
1442       return
1443       end
1444 c----------------------------------------------------------------------------
1445       subroutine csaread
1446       implicit real*8 (a-h,o-z)
1447       include 'DIMENSIONS'
1448       include 'COMMON.IOUNITS'
1449       include 'COMMON.GEO'
1450       include 'COMMON.CSA'
1451       include 'COMMON.BANK'
1452       include 'COMMON.CONTROL'
1453       character*80 ucase
1454       character*620 mcmcard
1455       call card_concat(mcmcard)
1456
1457       call readi(mcmcard,'NCONF',nconf,50)
1458       call readi(mcmcard,'NADD',nadd,0)
1459       call readi(mcmcard,'JSTART',jstart,1)
1460       call readi(mcmcard,'JEND',jend,1)
1461       call readi(mcmcard,'NSTMAX',nstmax,500000)
1462       call readi(mcmcard,'N0',n0,1)
1463       call readi(mcmcard,'N1',n1,6)
1464       call readi(mcmcard,'N2',n2,4)
1465       call readi(mcmcard,'N3',n3,0)
1466       call readi(mcmcard,'N4',n4,0)
1467       call readi(mcmcard,'N5',n5,0)
1468       call readi(mcmcard,'N6',n6,10)
1469       call readi(mcmcard,'N7',n7,0)
1470       call readi(mcmcard,'N8',n8,0)
1471       call readi(mcmcard,'N9',n9,0)
1472       call readi(mcmcard,'N14',n14,0)
1473       call readi(mcmcard,'N15',n15,0)
1474       call readi(mcmcard,'N16',n16,0)
1475       call readi(mcmcard,'N17',n17,0)
1476       call readi(mcmcard,'N18',n18,0)
1477
1478       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1479
1480       call readi(mcmcard,'NDIFF',ndiff,2)
1481       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1482       call readi(mcmcard,'IS1',is1,1)
1483       call readi(mcmcard,'IS2',is2,8)
1484       call readi(mcmcard,'NRAN0',nran0,4)
1485       call readi(mcmcard,'NRAN1',nran1,2)
1486       call readi(mcmcard,'IRR',irr,1)
1487       call readi(mcmcard,'NSEED',nseed,20)
1488       call readi(mcmcard,'NTOTAL',ntotal,10000)
1489       call reada(mcmcard,'CUT1',cut1,2.0d0)
1490       call reada(mcmcard,'CUT2',cut2,5.0d0)
1491       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1492       call readi(mcmcard,'ICMAX',icmax,3)
1493       call readi(mcmcard,'IRESTART',irestart,0)
1494 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1495       ntbankm=0
1496 c!bankt
1497       call reada(mcmcard,'DELE',dele,20.0d0)
1498       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1499       call readi(mcmcard,'IREF',iref,0)
1500       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1501       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1502       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1503       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1504       write (iout,*) "NCONF_IN",nconf_in
1505       return
1506       end
1507 c----------------------------------------------------------------------------
1508 cfmc      subroutine mcmfread
1509 cfmc      implicit real*8 (a-h,o-z)
1510 cfmc      include 'DIMENSIONS'
1511 cfmc      include 'COMMON.MCMF'
1512 cfmc      include 'COMMON.IOUNITS'
1513 cfmc      include 'COMMON.GEO'
1514 cfmc      character*80 ucase
1515 cfmc      character*620 mcmcard
1516 cfmc      call card_concat(mcmcard)
1517 cfmc
1518 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1519 cfmc      write(iout,*)'MAXRANT=',maxrant
1520 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1521 cfmc      write(iout,*)'MAXFAM=',maxfam
1522 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1523 cfmc      write(iout,*)'NNET1=',nnet1
1524 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1525 cfmc      write(iout,*)'NNET2=',nnet2
1526 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1527 cfmc      write(iout,*)'NNET3=',nnet3
1528 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1529 cfmc      write(iout,*)'ILASTT=',ilastt
1530 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1531 cfmc      write(iout,*)'MAXSTR=',maxstr
1532 cfmc      maxstr_f=maxstr/maxfam
1533 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1534 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1535 cfmc      write(iout,*)'NMCMF=',nmcmf
1536 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1537 cfmc      write(iout,*)'IFOCUS=',ifocus
1538 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1539 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1540 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1541 cfmc      write(iout,*)'INTPRT=',intprt
1542 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1543 cfmc      write(iout,*)'IPRT=',iprt
1544 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1545 cfmc      write(iout,*)'IMAXTR=',imaxtr
1546 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1547 cfmc      write(iout,*)'MAXEVEN=',maxeven
1548 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1549 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1550 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1551 cfmc      write(iout,*)'INIMIN=',inimin
1552 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1553 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1554 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1555 cfmc      write(iout,*)'NTHREAD=',nthread
1556 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1557 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1558 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1559 cfmc      write(iout,*)'MAXPERT=',maxpert
1560 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1561 cfmc      write(iout,*)'IRMSD=',irmsd
1562 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1563 cfmc      write(iout,*)'DENEMIN=',denemin
1564 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1565 cfmc      write(iout,*)'RCUT1S=',rcut1s
1566 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1567 cfmc      write(iout,*)'RCUT1E=',rcut1e
1568 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1569 cfmc      write(iout,*)'RCUT2S=',rcut2s
1570 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1571 cfmc      write(iout,*)'RCUT2E=',rcut2e
1572 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1573 cfmc      write(iout,*)'DPERT1=',d_pert1
1574 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1575 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1576 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1577 cfmc      write(iout,*)'DPERT2=',d_pert2
1578 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1579 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1580 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1581 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1582 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1583 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1584 cfmc      d_pert1=deg2rad*d_pert1
1585 cfmc      d_pert1a=deg2rad*d_pert1a
1586 cfmc      d_pert2=deg2rad*d_pert2
1587 cfmc      d_pert2a=deg2rad*d_pert2a
1588 cfmc      d_pert2b=deg2rad*d_pert2b
1589 cfmc      d_pert2c=deg2rad*d_pert2c
1590 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1591 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1592 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1593 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1594 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1595 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1596 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1597 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1598 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1599 cfmc      write(iout,*)'RCUTINI=',rcutini
1600 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1601 cfmc      write(iout,*)'GRAT=',grat
1602 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1603 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1604 cfmc
1605 cfmc      return
1606 cfmc      end
1607 c----------------------------------------------------------------------------
1608       subroutine mcmread
1609       implicit real*8 (a-h,o-z)
1610       include 'DIMENSIONS'
1611       include 'COMMON.MCM'
1612       include 'COMMON.MCE'
1613       include 'COMMON.IOUNITS'
1614       character*80 ucase
1615       character*320 mcmcard
1616       call card_concat(mcmcard)
1617       call readi(mcmcard,'MAXACC',maxacc,100)
1618       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1619       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1620       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1621       call readi(mcmcard,'MAXREPM',maxrepm,200)
1622       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1623       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1624       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1625       call reada(mcmcard,'E_UP',e_up,5.0D0)
1626       call reada(mcmcard,'DELTE',delte,0.1D0)
1627       call readi(mcmcard,'NSWEEP',nsweep,5)
1628       call readi(mcmcard,'NSTEPH',nsteph,0)
1629       call readi(mcmcard,'NSTEPC',nstepc,0)
1630       call reada(mcmcard,'TMIN',tmin,298.0D0)
1631       call reada(mcmcard,'TMAX',tmax,298.0D0)
1632       call readi(mcmcard,'NWINDOW',nwindow,0)
1633       call readi(mcmcard,'PRINT_MC',print_mc,0)
1634       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1635       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1636       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1637       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1638       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1639       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1640       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1641       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1642       if (nwindow.gt.0) then
1643         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1644         do i=1,nwindow
1645           winlen(i)=winend(i)-winstart(i)+1
1646         enddo
1647       endif
1648       if (tmax.lt.tmin) tmax=tmin
1649       if (tmax.eq.tmin) then
1650         nstepc=0
1651         nsteph=0
1652       endif
1653       if (nstepc.gt.0 .and. nsteph.gt.0) then
1654         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1655         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1656       endif
1657 C Probabilities of different move types
1658       sumpro_type(0)=0.0D0
1659       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1660       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1661       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1662       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1663       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1664       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1665       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1666       do i=1,MaxMoveType
1667         print *,'i',i,' sumprotype',sumpro_type(i)
1668         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1669         print *,'i',i,' sumprotype',sumpro_type(i)
1670       enddo
1671       return
1672       end
1673 c----------------------------------------------------------------------------
1674       subroutine read_minim
1675       implicit real*8 (a-h,o-z)
1676       include 'DIMENSIONS'
1677       include 'COMMON.MINIM'
1678       include 'COMMON.IOUNITS'
1679       character*80 ucase
1680       character*320 minimcard
1681       call card_concat(minimcard)
1682       call readi(minimcard,'MAXMIN',maxmin,2000)
1683       call readi(minimcard,'MAXFUN',maxfun,5000)
1684       call readi(minimcard,'MINMIN',minmin,maxmin)
1685       call readi(minimcard,'MINFUN',minfun,maxmin)
1686       call reada(minimcard,'TOLF',tolf,1.0D-2)
1687       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1688       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1689       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1690       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1691       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1692      &         'Options in energy minimization:'
1693       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1694      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1695      & 'MinMin:',MinMin,' MinFun:',MinFun,
1696      & ' TolF:',TolF,' RTolF:',RTolF
1697       return
1698       end
1699 c----------------------------------------------------------------------------
1700       subroutine read_angles(kanal,*)
1701       implicit real*8 (a-h,o-z)
1702       include 'DIMENSIONS'
1703       include 'COMMON.GEO'
1704       include 'COMMON.VAR'
1705       include 'COMMON.CHAIN'
1706       include 'COMMON.IOUNITS'
1707       include 'COMMON.CONTROL'
1708 c Read angles from input
1709 c
1710        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1711        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1712        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1713        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1714
1715        do i=1,nres
1716 c 9/7/01 avoid 180 deg valence angle
1717         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1718 c
1719         theta(i)=deg2rad*theta(i)
1720         phi(i)=deg2rad*phi(i)
1721         alph(i)=deg2rad*alph(i)
1722         omeg(i)=deg2rad*omeg(i)
1723        enddo
1724       return
1725    10 return1
1726       end
1727 c----------------------------------------------------------------------------
1728       subroutine reada(rekord,lancuch,wartosc,default)
1729       implicit none
1730       character*(*) rekord,lancuch
1731       double precision wartosc,default
1732       integer ilen,iread
1733       external ilen
1734       iread=index(rekord,lancuch)
1735       if (iread.eq.0) then
1736         wartosc=default
1737         return
1738       endif
1739       iread=iread+ilen(lancuch)+1
1740       read (rekord(iread:),*,err=10,end=10) wartosc
1741       return
1742   10  wartosc=default
1743       return
1744       end
1745 c----------------------------------------------------------------------------
1746       subroutine readi(rekord,lancuch,wartosc,default)
1747       implicit none
1748       character*(*) rekord,lancuch
1749       integer wartosc,default
1750       integer ilen,iread
1751       external ilen
1752       iread=index(rekord,lancuch)
1753       if (iread.eq.0) then
1754         wartosc=default
1755         return
1756       endif
1757       iread=iread+ilen(lancuch)+1
1758       read (rekord(iread:),*,err=10,end=10) wartosc
1759       return
1760   10  wartosc=default
1761       return
1762       end
1763 c----------------------------------------------------------------------------
1764       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1765       implicit none
1766       integer dim,i
1767       integer tablica(dim),default
1768       character*(*) rekord,lancuch
1769       character*80 aux
1770       integer ilen,iread
1771       external ilen
1772       do i=1,dim
1773         tablica(i)=default
1774       enddo
1775       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1776       if (iread.eq.0) return
1777       iread=iread+ilen(lancuch)+1
1778       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1779    10 return
1780       end
1781 c----------------------------------------------------------------------------
1782       subroutine multreada(rekord,lancuch,tablica,dim,default)
1783       implicit none
1784       integer dim,i
1785       double precision tablica(dim),default
1786       character*(*) rekord,lancuch
1787       character*80 aux
1788       integer ilen,iread
1789       external ilen
1790       do i=1,dim
1791         tablica(i)=default
1792       enddo
1793       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1794       if (iread.eq.0) return
1795       iread=iread+ilen(lancuch)+1
1796       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1797    10 return
1798       end
1799 c----------------------------------------------------------------------------
1800       subroutine openunits
1801       implicit real*8 (a-h,o-z)
1802       include 'DIMENSIONS'
1803 #ifdef MPI
1804       include 'mpif.h'
1805       character*16 form,nodename
1806       integer nodelen
1807 #endif
1808       include 'COMMON.SETUP'
1809       include 'COMMON.IOUNITS'
1810       include 'COMMON.MD'
1811       include 'COMMON.CONTROL'
1812       integer lenpre,lenpot,ilen,lentmp
1813       external ilen
1814       character*3 out1file_text,ucase
1815       character*3 ll
1816       external ucase
1817 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1818       call getenv_loc("PREFIX",prefix)
1819       pref_orig = prefix
1820       call getenv_loc("POT",pot)
1821       call getenv_loc("DIRTMP",tmpdir)
1822       call getenv_loc("CURDIR",curdir)
1823       call getenv_loc("OUT1FILE",out1file_text)
1824 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1825       out1file_text=ucase(out1file_text)
1826       if (out1file_text(1:1).eq."Y") then
1827         out1file=.true.
1828       else
1829         out1file=fg_rank.gt.0
1830       endif
1831       lenpre=ilen(prefix)
1832       lenpot=ilen(pot)
1833       lentmp=ilen(tmpdir)
1834       if (lentmp.gt.0) then
1835           write (*,'(80(1h!))')
1836           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1837           write (*,'(80(1h!))')
1838           write (*,*)"All output files will be on node /tmp directory."
1839 #ifdef MPI
1840         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1841         if (me.eq.king) then
1842           write (*,*) "The master node is ",nodename
1843         else if (fg_rank.eq.0) then
1844           write (*,*) "I am the CG slave node ",nodename
1845         else
1846           write (*,*) "I am the FG slave node ",nodename
1847         endif
1848 #endif
1849         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1850         lenpre = lentmp+lenpre+1
1851       endif
1852       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1853 C Get the names and open the input files
1854 #if defined(WINIFL) || defined(WINPGI)
1855       open(1,file=pref_orig(:ilen(pref_orig))//
1856      &  '.inp',status='old',readonly,shared)
1857        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1858 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1859 C Get parameter filenames and open the parameter files.
1860       call getenv_loc('BONDPAR',bondname)
1861       open (ibond,file=bondname,status='old',readonly,shared)
1862       call getenv_loc('THETPAR',thetname)
1863       open (ithep,file=thetname,status='old',readonly,shared)
1864       call getenv_loc('ROTPAR',rotname)
1865       open (irotam,file=rotname,status='old',readonly,shared)
1866       call getenv_loc('TORPAR',torname)
1867       open (itorp,file=torname,status='old',readonly,shared)
1868       call getenv_loc('TORDPAR',tordname)
1869       open (itordp,file=tordname,status='old',readonly,shared)
1870       call getenv_loc('FOURIER',fouriername)
1871       open (ifourier,file=fouriername,status='old',readonly,shared)
1872       call getenv_loc('ELEPAR',elename)
1873       open (ielep,file=elename,status='old',readonly,shared)
1874       call getenv_loc('SIDEPAR',sidename)
1875       open (isidep,file=sidename,status='old',readonly,shared)
1876 #elif (defined CRAY) || (defined AIX)
1877       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1878      &  action='read')
1879 c      print *,"Processor",myrank," opened file 1"
1880       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1881 c      print *,"Processor",myrank," opened file 9"
1882 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1883 C Get parameter filenames and open the parameter files.
1884       call getenv_loc('BONDPAR',bondname)
1885       open (ibond,file=bondname,status='old',action='read')
1886 c      print *,"Processor",myrank," opened file IBOND"
1887       call getenv_loc('THETPAR',thetname)
1888       open (ithep,file=thetname,status='old',action='read')
1889 c      print *,"Processor",myrank," opened file ITHEP"
1890       call getenv_loc('ROTPAR',rotname)
1891       open (irotam,file=rotname,status='old',action='read')
1892 c      print *,"Processor",myrank," opened file IROTAM"
1893       call getenv_loc('TORPAR',torname)
1894       open (itorp,file=torname,status='old',action='read')
1895 c      print *,"Processor",myrank," opened file ITORP"
1896       call getenv_loc('TORDPAR',tordname)
1897       open (itordp,file=tordname,status='old',action='read')
1898 c      print *,"Processor",myrank," opened file ITORDP"
1899       call getenv_loc('SCCORPAR',sccorname)
1900       open (isccor,file=sccorname,status='old',action='read')
1901 c      print *,"Processor",myrank," opened file ISCCOR"
1902       call getenv_loc('FOURIER',fouriername)
1903       open (ifourier,file=fouriername,status='old',action='read')
1904 c      print *,"Processor",myrank," opened file IFOURIER"
1905       call getenv_loc('ELEPAR',elename)
1906       open (ielep,file=elename,status='old',action='read')
1907 c      print *,"Processor",myrank," opened file IELEP"
1908       call getenv_loc('SIDEPAR',sidename)
1909       open (isidep,file=sidename,status='old',action='read')
1910 c      print *,"Processor",myrank," opened file ISIDEP"
1911 c      print *,"Processor",myrank," opened parameter files"
1912 #elif (defined G77)
1913       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1914       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1915 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1916 C Get parameter filenames and open the parameter files.
1917       call getenv_loc('BONDPAR',bondname)
1918       open (ibond,file=bondname,status='old')
1919       call getenv_loc('THETPAR',thetname)
1920       open (ithep,file=thetname,status='old')
1921       call getenv_loc('ROTPAR',rotname)
1922       open (irotam,file=rotname,status='old')
1923       call getenv_loc('TORPAR',torname)
1924       open (itorp,file=torname,status='old')
1925       call getenv_loc('TORDPAR',tordname)
1926       open (itordp,file=tordname,status='old')
1927       call getenv_loc('SCCORPAR',sccorname)
1928       open (isccor,file=sccorname,status='old')
1929       call getenv_loc('FOURIER',fouriername)
1930       open (ifourier,file=fouriername,status='old')
1931       call getenv_loc('ELEPAR',elename)
1932       open (ielep,file=elename,status='old')
1933       call getenv_loc('SIDEPAR',sidename)
1934       open (isidep,file=sidename,status='old')
1935 #else
1936       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1937      &  readonly)
1938        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1939 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1940 C Get parameter filenames and open the parameter files.
1941       call getenv_loc('BONDPAR',bondname)
1942       open (ibond,file=bondname,status='old',readonly)
1943       call getenv_loc('THETPAR',thetname)
1944       open (ithep,file=thetname,status='old',readonly)
1945       call getenv_loc('ROTPAR',rotname)
1946       open (irotam,file=rotname,status='old',readonly)
1947       call getenv_loc('TORPAR',torname)
1948       open (itorp,file=torname,status='old',readonly)
1949       call getenv_loc('TORDPAR',tordname)
1950       open (itordp,file=tordname,status='old',readonly)
1951       call getenv_loc('SCCORPAR',sccorname)
1952       open (isccor,file=sccorname,status='old',readonly)
1953 #ifndef CRYST_THETA
1954       call getenv_loc('THETPARPDB',thetname_pdb)
1955       print *,"thetname_pdb ",thetname_pdb
1956       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1957       print *,ithep_pdb," opened"
1958 #endif
1959       call getenv_loc('FOURIER',fouriername)
1960       open (ifourier,file=fouriername,status='old',readonly)
1961       call getenv_loc('ELEPAR',elename)
1962       open (ielep,file=elename,status='old',readonly)
1963       call getenv_loc('SIDEPAR',sidename)
1964       open (isidep,file=sidename,status='old',readonly)
1965 #ifndef CRYST_SC
1966       call getenv_loc('ROTPARPDB',rotname_pdb)
1967       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1968 #endif
1969 #endif
1970 #ifndef OLDSCP
1971 C
1972 C 8/9/01 In the newest version SCp interaction constants are read from a file
1973 C Use -DOLDSCP to use hard-coded constants instead.
1974 C
1975       call getenv_loc('SCPPAR',scpname)
1976 #if defined(WINIFL) || defined(WINPGI)
1977       open (iscpp,file=scpname,status='old',readonly,shared)
1978 #elif (defined CRAY)  || (defined AIX)
1979       open (iscpp,file=scpname,status='old',action='read')
1980 #elif (defined G77)
1981       open (iscpp,file=scpname,status='old')
1982 #else
1983       open (iscpp,file=scpname,status='old',readonly)
1984 #endif
1985 #endif
1986       call getenv_loc('PATTERN',patname)
1987 #if defined(WINIFL) || defined(WINPGI)
1988       open (icbase,file=patname,status='old',readonly,shared)
1989 #elif (defined CRAY)  || (defined AIX)
1990       open (icbase,file=patname,status='old',action='read')
1991 #elif (defined G77)
1992       open (icbase,file=patname,status='old')
1993 #else
1994       open (icbase,file=patname,status='old',readonly)
1995 #endif
1996 #ifdef MPI
1997 C Open output file only for CG processes
1998 c      print *,"Processor",myrank," fg_rank",fg_rank
1999       if (fg_rank.eq.0) then
2000
2001       if (nodes.eq.1) then
2002         npos=3
2003       else
2004         npos = dlog10(dfloat(nodes-1))+1
2005       endif
2006       if (npos.lt.3) npos=3
2007       write (liczba,'(i1)') npos
2008       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2009      &  //')'
2010       write (liczba,form) me
2011       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2012      &  liczba(:ilen(liczba))
2013       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2014      &  //'.int'
2015       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2016      &  //'.pdb'
2017       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2018      &  liczba(:ilen(liczba))//'.mol2'
2019       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2020      &  liczba(:ilen(liczba))//'.stat'
2021       if (lentmp.gt.0)
2022      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2023      &      //liczba(:ilen(liczba))//'.stat')
2024       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2025      &  //'.rst'
2026       if(usampl) then
2027           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2028      & liczba(:ilen(liczba))//'.const'
2029       endif
2030
2031       endif
2032 #else
2033       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2034       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2035       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2036       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2037       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2038       if (lentmp.gt.0)
2039      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2040      &    //'.stat')
2041       rest2name=prefix(:ilen(prefix))//'.rst'
2042       if(usampl) then
2043          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2044       endif
2045 #endif
2046 #if defined(AIX) || defined(PGI)
2047       if (me.eq.king .or. .not. out1file)
2048      &   open(iout,file=outname,status='unknown')
2049 #ifdef DEBUG
2050       if (fg_rank.gt.0) then
2051         write (liczba,'(i3.3)') myrank/nfgtasks
2052         write (ll,'(bz,i3.3)') fg_rank
2053         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2054      &   status='unknown')
2055       endif
2056 #endif
2057       if(me.eq.king) then
2058        open(igeom,file=intname,status='unknown',position='append')
2059        open(ipdb,file=pdbname,status='unknown')
2060        open(imol2,file=mol2name,status='unknown')
2061        open(istat,file=statname,status='unknown',position='append')
2062       else
2063 c1out       open(iout,file=outname,status='unknown')
2064       endif
2065 #else
2066       if (me.eq.king .or. .not.out1file)
2067      &    open(iout,file=outname,status='unknown')
2068 #ifdef DEBUG
2069       if (fg_rank.gt.0) then
2070         write (liczba,'(i3.3)') myrank/nfgtasks
2071         write (ll,'(bz,i3.3)') fg_rank
2072         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2073      &   status='unknown')
2074       endif
2075 #endif
2076       if(me.eq.king) then
2077        open(igeom,file=intname,status='unknown',access='append')
2078        open(ipdb,file=pdbname,status='unknown')
2079        open(imol2,file=mol2name,status='unknown')
2080        open(istat,file=statname,status='unknown',access='append')
2081       else
2082 c1out       open(iout,file=outname,status='unknown')
2083       endif
2084 #endif
2085       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2086       csa_seed=prefix(:lenpre)//'.CSA.seed'
2087       csa_history=prefix(:lenpre)//'.CSA.history'
2088       csa_bank=prefix(:lenpre)//'.CSA.bank'
2089       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2090       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2091       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2092 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2093       csa_int=prefix(:lenpre)//'.int'
2094       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2095       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2096       csa_in=prefix(:lenpre)//'.CSA.in'
2097 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2098 C Write file names
2099       if (me.eq.king)then
2100       write (iout,'(80(1h-))')
2101       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2102       write (iout,'(80(1h-))')
2103       write (iout,*) "Input file                      : ",
2104      &  pref_orig(:ilen(pref_orig))//'.inp'
2105       write (iout,*) "Output file                     : ",
2106      &  outname(:ilen(outname))
2107       write (iout,*)
2108       write (iout,*) "Sidechain potential file        : ",
2109      &  sidename(:ilen(sidename))
2110 #ifndef OLDSCP
2111       write (iout,*) "SCp potential file              : ",
2112      &  scpname(:ilen(scpname))
2113 #endif
2114       write (iout,*) "Electrostatic potential file    : ",
2115      &  elename(:ilen(elename))
2116       write (iout,*) "Cumulant coefficient file       : ",
2117      &  fouriername(:ilen(fouriername))
2118       write (iout,*) "Torsional parameter file        : ",
2119      &  torname(:ilen(torname))
2120       write (iout,*) "Double torsional parameter file : ",
2121      &  tordname(:ilen(tordname))
2122       write (iout,*) "SCCOR parameter file : ",
2123      &  sccorname(:ilen(sccorname))
2124       write (iout,*) "Bond & inertia constant file    : ",
2125      &  bondname(:ilen(bondname))
2126       write (iout,*) "Bending parameter file          : ",
2127      &  thetname(:ilen(thetname))
2128       write (iout,*) "Rotamer parameter file          : ",
2129      &  rotname(:ilen(rotname))
2130       write (iout,*) "Thetpdb parameter file          : ",
2131      &  thetname_pdb(:ilen(thetname_pdb))
2132       write (iout,*) "Threading database              : ",
2133      &  patname(:ilen(patname))
2134       if (lentmp.ne.0)
2135      &write (iout,*)" DIRTMP                          : ",
2136      &  tmpdir(:lentmp)
2137       write (iout,'(80(1h-))')
2138       endif
2139       return
2140       end
2141 c----------------------------------------------------------------------------
2142       subroutine card_concat(card)
2143       implicit real*8 (a-h,o-z)
2144       include 'DIMENSIONS'
2145       include 'COMMON.IOUNITS'
2146       character*(*) card
2147       character*80 karta,ucase
2148       external ilen
2149       read (inp,'(a)') karta
2150       karta=ucase(karta)
2151       card=' '
2152       do while (karta(80:80).eq.'&')
2153         card=card(:ilen(card)+1)//karta(:79)
2154         read (inp,'(a)') karta
2155         karta=ucase(karta)
2156       enddo
2157       card=card(:ilen(card)+1)//karta
2158       return
2159       end
2160 c----------------------------------------------------------------------------------
2161       subroutine readrst
2162       implicit real*8 (a-h,o-z)
2163       include 'DIMENSIONS'
2164       include 'COMMON.CHAIN'
2165       include 'COMMON.IOUNITS'
2166       include 'COMMON.MD'
2167       open(irest2,file=rest2name,status='unknown')
2168       read(irest2,*) totT,EK,potE,totE,t_bath
2169       do i=1,2*nres
2170          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2171       enddo
2172       do i=1,2*nres
2173          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2174       enddo
2175       if(usampl) then
2176              read (irest2,*) iset
2177       endif
2178       close(irest2)
2179       return
2180       end
2181 c---------------------------------------------------------------------------------
2182       subroutine read_fragments
2183       implicit real*8 (a-h,o-z)
2184       include 'DIMENSIONS'
2185 #ifdef MPI
2186       include 'mpif.h'
2187 #endif
2188       include 'COMMON.SETUP'
2189       include 'COMMON.CHAIN'
2190       include 'COMMON.IOUNITS'
2191       include 'COMMON.MD'
2192       include 'COMMON.CONTROL'
2193       read(inp,*) nset,nfrag,npair,nfrag_back
2194       if(me.eq.king.or..not.out1file)
2195      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2196      &  " nfrag_back",nfrag_back
2197       do iset=1,nset
2198          read(inp,*) mset(iset)
2199        do i=1,nfrag
2200          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2201      &     qinfrag(i,iset)
2202          if(me.eq.king.or..not.out1file)
2203      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2204      &     ifrag(2,i,iset), qinfrag(i,iset)
2205        enddo
2206        do i=1,npair
2207         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2208      &    qinpair(i,iset)
2209         if(me.eq.king.or..not.out1file)
2210      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2211      &    ipair(2,i,iset), qinpair(i,iset)
2212        enddo
2213        do i=1,nfrag_back
2214         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2215      &     wfrag_back(3,i,iset),
2216      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2217         if(me.eq.king.or..not.out1file)
2218      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2219      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2220        enddo
2221       enddo
2222       return
2223       end
2224 c-------------------------------------------------------------------------------
2225       subroutine read_dist_constr
2226       implicit real*8 (a-h,o-z)
2227       include 'DIMENSIONS'
2228 #ifdef MPI
2229       include 'mpif.h'
2230 #endif
2231       include 'COMMON.SETUP'
2232       include 'COMMON.CONTROL'
2233       include 'COMMON.CHAIN'
2234       include 'COMMON.IOUNITS'
2235       include 'COMMON.SBRIDGE'
2236       integer ifrag_(2,100),ipair_(2,100)
2237       double precision wfrag_(100),wpair_(100)
2238       character*500 controlcard
2239 c      write (iout,*) "Calling read_dist_constr"
2240 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2241 c      call flush(iout)
2242       call card_concat(controlcard)
2243       call readi(controlcard,"NFRAG",nfrag_,0)
2244       call readi(controlcard,"NPAIR",npair_,0)
2245       call readi(controlcard,"NDIST",ndist_,0)
2246       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2247       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2248       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2249       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2250       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2251 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2252 c      write (iout,*) "IFRAG"
2253 c      do i=1,nfrag_
2254 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2255 c      enddo
2256 c      write (iout,*) "IPAIR"
2257 c      do i=1,npair_
2258 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2259 c      enddo
2260       call flush(iout)
2261       do i=1,nfrag_
2262         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2263         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2264      &    ifrag_(2,i)=nstart_sup+nsup-1
2265 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2266         call flush(iout)
2267         if (wfrag_(i).gt.0.0d0) then
2268         do j=ifrag_(1,i),ifrag_(2,i)-1
2269           do k=j+1,ifrag_(2,i)
2270 c            write (iout,*) "j",j," k",k
2271             ddjk=dist(j,k)
2272             if (constr_dist.eq.1) then
2273             nhpb=nhpb+1
2274             ihpb(nhpb)=j
2275             jhpb(nhpb)=k
2276               dhpb(nhpb)=ddjk
2277             forcon(nhpb)=wfrag_(i)
2278             else if (constr_dist.eq.2) then
2279               if (ddjk.le.dist_cut) then
2280                 nhpb=nhpb+1
2281                 ihpb(nhpb)=j
2282                 jhpb(nhpb)=k
2283                 dhpb(nhpb)=ddjk
2284                 forcon(nhpb)=wfrag_(i)
2285               endif
2286             else
2287               nhpb=nhpb+1
2288               ihpb(nhpb)=j
2289               jhpb(nhpb)=k
2290               dhpb(nhpb)=ddjk
2291               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2292             endif
2293 #ifdef MPI
2294             if (.not.out1file .or. me.eq.king)
2295      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2296      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2297 #else
2298             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2299      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2300 #endif
2301           enddo
2302         enddo
2303         endif
2304       enddo
2305       do i=1,npair_
2306         if (wpair_(i).gt.0.0d0) then
2307         ii = ipair_(1,i)
2308         jj = ipair_(2,i)
2309         if (ii.gt.jj) then
2310           itemp=ii
2311           ii=jj
2312           jj=itemp
2313         endif
2314         do j=ifrag_(1,ii),ifrag_(2,ii)
2315           do k=ifrag_(1,jj),ifrag_(2,jj)
2316             nhpb=nhpb+1
2317             ihpb(nhpb)=j
2318             jhpb(nhpb)=k
2319             forcon(nhpb)=wpair_(i)
2320             dhpb(nhpb)=dist(j,k)
2321 #ifdef MPI
2322             if (.not.out1file .or. me.eq.king)
2323      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2324      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2325 #else
2326             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2327      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2328 #endif
2329           enddo
2330         enddo
2331         endif
2332       enddo
2333       do i=1,ndist_
2334         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2335         if (forcon(nhpb+1).gt.0.0d0) then
2336           nhpb=nhpb+1
2337           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2338 #ifdef MPI
2339           if (.not.out1file .or. me.eq.king)
2340      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2341      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2342 #else
2343           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2344      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2345 #endif
2346         endif
2347       enddo
2348       call flush(iout)
2349       return
2350       end
2351 c-------------------------------------------------------------------------------
2352 #ifdef WINIFL
2353       subroutine flush(iu)
2354       return
2355       end
2356 #endif
2357 #ifdef AIX
2358       subroutine flush(iu)
2359       call flush_(iu)
2360       return
2361       end
2362 #endif
2363 c------------------------------------------------------------------------------
2364       subroutine copy_to_tmp(source)
2365       include "DIMENSIONS"
2366       include "COMMON.IOUNITS"
2367       character*(*) source
2368       character* 256 tmpfile
2369       integer ilen
2370       external ilen
2371       logical ex
2372       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2373       inquire(file=tmpfile,exist=ex)
2374       if (ex) then
2375         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2376      &   " to temporary directory..."
2377         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2378         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2379       endif
2380       return
2381       end
2382 c------------------------------------------------------------------------------
2383       subroutine move_from_tmp(source)
2384       include "DIMENSIONS"
2385       include "COMMON.IOUNITS"
2386       character*(*) source
2387       integer ilen
2388       external ilen
2389       write (*,*) "Moving ",source(:ilen(source)),
2390      & " from temporary directory to working directory"
2391       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2392       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2393       return
2394       end
2395 c------------------------------------------------------------------------------
2396       subroutine random_init(seed)
2397 C
2398 C Initialize random number generator
2399 C
2400       implicit real*8 (a-h,o-z)
2401       include 'DIMENSIONS'
2402 #ifdef MPI
2403       include 'mpif.h'
2404       logical OKRandom, prng_restart
2405       real*8  r1
2406       integer iseed_array(4)
2407 #endif
2408       include 'COMMON.IOUNITS'
2409       include 'COMMON.TIME1'
2410       include 'COMMON.THREAD'
2411       include 'COMMON.SBRIDGE'
2412       include 'COMMON.CONTROL'
2413       include 'COMMON.MCM'
2414       include 'COMMON.MAP'
2415       include 'COMMON.HEADER'
2416       include 'COMMON.CSA'
2417       include 'COMMON.CHAIN'
2418       include 'COMMON.MUCA'
2419       include 'COMMON.MD'
2420       include 'COMMON.FFIELD'
2421       include 'COMMON.SETUP'
2422       iseed=-dint(dabs(seed))
2423       if (iseed.eq.0) then
2424         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2425      &    'Random seed undefined. The program will stop.'
2426         write (*,'(/80(1h*)/20x,a/80(1h*))')
2427      &    'Random seed undefined. The program will stop.'
2428 #ifdef MPI
2429         call mpi_finalize(mpi_comm_world,ierr)
2430 #endif
2431         stop 'Bad random seed.'
2432       endif
2433 #ifdef MPI
2434       if (fg_rank.eq.0) then
2435       seed=seed*(me+1)+1
2436 #ifdef AMD64
2437       if(me.eq.king)
2438      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2439       OKRandom = prng_restart(me,iseed)
2440 #else
2441       do i=1,4
2442        tmp=65536.0d0**(4-i)
2443        iseed_array(i) = dint(seed/tmp)
2444        seed=seed-iseed_array(i)*tmp
2445       enddo
2446       if(me.eq.king)
2447      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2448      &                 (iseed_array(i),i=1,4)
2449       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2450      &                 (iseed_array(i),i=1,4)
2451       OKRandom = prng_restart(me,iseed_array)
2452 #endif
2453       if (OKRandom) then
2454 c        r1 = prng_next(me)
2455         r1=ran_number(0.0D0,1.0D0)
2456         if(me.eq.king)
2457      &   write (iout,*) 'ran_num',r1
2458         if (r1.lt.0.0d0) OKRandom=.false.
2459       endif
2460       if (.not.OKRandom) then
2461         write (iout,*) 'PRNG IS NOT WORKING!!!'
2462         print *,'PRNG IS NOT WORKING!!!'
2463         if (me.eq.0) then
2464          call flush(iout)
2465          call mpi_abort(mpi_comm_world,error_msg,ierr)
2466          stop
2467         else
2468          write (iout,*) 'too many processors for parallel prng'
2469          write (*,*) 'too many processors for parallel prng'
2470          call flush(iout)
2471          stop
2472         endif
2473       endif
2474       endif
2475 #else
2476       call vrndst(iseed)
2477       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2478 #endif
2479       return
2480       end