--- /dev/null
+ subroutine readrtns
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.IOUNITS'
+ logical file_exist
+C Read force-field parameters except weights
+ call parmread
+C Read job setup parameters
+ call read_control
+C Read control parameters for energy minimzation if required
+ if (minim) call read_minim
+C Read MCM control parameters if required
+ if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
+C Read MD control parameters if reqjuired
+ if (modecalc.eq.12) call read_MDpar
+C Read MREMD control parameters if required
+ if (modecalc.eq.14) then
+ call read_MDpar
+ call read_REMDpar
+ endif
+C Read MUCA control parameters if required
+ if (lmuca) call read_muca
+C Read CSA control parameters if required (from fort.40 if exists
+C otherwise from general input file)
+csa if (modecalc.eq.8) then
+csa inquire (file="fort.40",exist=file_exist)
+csa if (.not.file_exist) call csaread
+csa endif
+cfmc if (modecalc.eq.10) call mcmfread
+C Read molecule information, molecule geometry, energy-term weights, and
+C restraints if requested
+ call molread
+C Print restraint information
+#ifdef MPI
+ if (.not. out1file .or. me.eq.king) then
+#endif
+ if (nhpb.gt.nss)
+ &write (iout,'(a,i5,a)') "The following",nhpb-nss,
+ & " distance constraints have been imposed"
+ do i=nss+1,nhpb
+ write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+c print *,"Processor",myrank," leaves READRTNS"
+ return
+ end
+C-------------------------------------------------------------------------------
+ subroutine read_control
+C
+C Read contorl data
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MP
+ include 'mpif.h'
+ logical OKRandom, prng_restart
+ real*8 r1
+#endif
+ include 'COMMON.IOUNITS'
+ include 'COMMON.TIME1'
+ include 'COMMON.THREAD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CONTROL'
+ include 'COMMON.MCM'
+ include 'COMMON.MAP'
+ include 'COMMON.HEADER'
+csa include 'COMMON.CSA'
+ include 'COMMON.CHAIN'
+ include 'COMMON.MUCA'
+ include 'COMMON.MD'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SETUP'
+ COMMON /MACHSW/ KDIAG,ICORFL,IXDR
+ character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
+ character*80 ucase
+ character*320 controlcard
+
+ nglob_csa=0
+ eglob_csa=1d99
+ nmin_csa=0
+ read (INP,'(a)') titel
+ call card_concat(controlcard)
+c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
+c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
+ call reada(controlcard,'SEED',seed,0.0D0)
+ call random_init(seed)
+C Set up the time limit (caution! The time must be input in minutes!)
+ read_cart=index(controlcard,'READ_CART').gt.0
+ call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
+ unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
+ call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
+ call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
+ call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
+ call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
+ call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
+ call reada(controlcard,'DRMS',drms,0.1D0)
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
+ write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
+ write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
+ write (iout,'(a,f10.1)')'DRMS = ',drms
+ write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
+ write (iout,'(a,f10.1)') 'Time limit (min):',timlim
+ endif
+ call readi(controlcard,'NZ_START',nz_start,0)
+ call readi(controlcard,'NZ_END',nz_end,0)
+ call readi(controlcard,'IZ_SC',iz_sc,0)
+ timlim=60.0D0*timlim
+ safety = 60.0d0*safety
+ timem=timlim
+ modecalc=0
+ call reada(controlcard,"T_BATH",t_bath,300.0d0)
+ minim=(index(controlcard,'MINIMIZE').gt.0)
+ dccart=(index(controlcard,'CART').gt.0)
+ overlapsc=(index(controlcard,'OVERLAP').gt.0)
+ overlapsc=.not.overlapsc
+ searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
+ searchsc=.not.searchsc
+ sideadd=(index(controlcard,'SIDEADD').gt.0)
+ energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
+ outpdb=(index(controlcard,'PDBOUT').gt.0)
+ outmol2=(index(controlcard,'MOL2OUT').gt.0)
+ pdbref=(index(controlcard,'PDBREF').gt.0)
+ refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
+ indpdb=index(controlcard,'PDBSTART')
+ extconf=(index(controlcard,'EXTCONF').gt.0)
+ call readi(controlcard,'IPRINT',iprint,0)
+ call readi(controlcard,'MAXGEN',maxgen,10000)
+ call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
+ call readi(controlcard,"KDIAG",kdiag,0)
+ call readi(controlcard,"RESCALE_MODE",rescale_mode,1)
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
+ & write (iout,*) "RESCALE_MODE",rescale_mode
+ split_ene=index(controlcard,'SPLIT_ENE').gt.0
+ if (index(controlcard,'REGULAR').gt.0.0D0) then
+ call reada(controlcard,'WEIDIS',weidis,0.1D0)
+ modecalc=1
+ refstr=.true.
+ endif
+ if (index(controlcard,'CHECKGRAD').gt.0) then
+ modecalc=5
+ if (index(controlcard,'CART').gt.0) then
+ icheckgrad=1
+ elseif (index(controlcard,'CARINT').gt.0) then
+ icheckgrad=2
+ else
+ icheckgrad=3
+ endif
+ elseif (index(controlcard,'THREAD').gt.0) then
+ modecalc=2
+ call readi(controlcard,'THREAD',nthread,0)
+ if (nthread.gt.0) then
+ call reada(controlcard,'WEIDIS',weidis,0.1D0)
+ else
+ if (fg_rank.eq.0)
+ & write (iout,'(a)')'A number has to follow the THREAD keyword.'
+ stop 'Error termination in Read_Control.'
+ endif
+ else if (index(controlcard,'MCMA').gt.0) then
+ modecalc=3
+ else if (index(controlcard,'MCEE').gt.0) then
+ modecalc=6
+ else if (index(controlcard,'MULTCONF').gt.0) then
+ modecalc=4
+ else if (index(controlcard,'MAP').gt.0) then
+ modecalc=7
+ call readi(controlcard,'MAP',nmap,0)
+ else if (index(controlcard,'CSA').gt.0) then
+ write(*,*) "CSA not supported in this version"
+ stop
+csa modecalc=8
+crc else if (index(controlcard,'ZSCORE').gt.0) then
+crc
+crc ZSCORE is rm from UNRES, modecalc=9 is available
+crc
+crc modecalc=9
+cfcm else if (index(controlcard,'MCMF').gt.0) then
+cfmc modecalc=10
+ else if (index(controlcard,'SOFTREG').gt.0) then
+ modecalc=11
+ else if (index(controlcard,'CHECK_BOND').gt.0) then
+ modecalc=-1
+ else if (index(controlcard,'TEST').gt.0) then
+ modecalc=-2
+ else if (index(controlcard,'MD').gt.0) then
+ modecalc=12
+ else if (index(controlcard,'RE ').gt.0) then
+ modecalc=14
+ endif
+
+ lmuca=index(controlcard,'MUCA').gt.0
+ call readi(controlcard,'MUCADYN',mucadyn,0)
+ call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
+ if (lmuca .and. (me.eq.king .or. .not.out1file ))
+ & then
+ write (iout,*) 'MUCADYN=',mucadyn
+ write (iout,*) 'MUCASMOOTH=',muca_smooth
+ endif
+
+ iscode=index(controlcard,'ONE_LETTER')
+ indphi=index(controlcard,'PHI')
+ indback=index(controlcard,'BACK')
+ iranconf=index(controlcard,'RAND_CONF')
+ i2ndstr=index(controlcard,'USE_SEC_PRED')
+ gradout=index(controlcard,'GRADOUT').gt.0
+ gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
+
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(2a)') diagmeth(kdiag),
+ & ' routine used to diagonalize matrices.'
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine read_REMDpar
+C
+C Read REMD settings
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.TIME1'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
+ include 'COMMON.GEO'
+ include 'COMMON.REMD'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SETUP'
+ character*80 ucase
+ character*320 controlcard
+ character*3200 controlcard1
+ integer iremd_m_total
+
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) "REMD setup"
+
+ call card_concat(controlcard)
+ call readi(controlcard,"NREP",nrep,3)
+ call readi(controlcard,"NSTEX",nstex,1000)
+ call reada(controlcard,"RETMIN",retmin,10.0d0)
+ call reada(controlcard,"RETMAX",retmax,1000.0d0)
+ mremdsync=(index(controlcard,'SYNC').gt.0)
+ call readi(controlcard,"NSYN",i_sync_step,100)
+ restart1file=(index(controlcard,'REST1FILE').gt.0)
+ traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
+ call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
+ if(max_cache_traj_use.gt.max_cache_traj)
+ & max_cache_traj_use=max_cache_traj
+ if(me.eq.king.or..not.out1file) then
+cd if (traj1file) then
+crc caching is in testing - NTWX is not ignored
+cd write (iout,*) "NTWX value is ignored"
+cd write (iout,*) " trajectory is stored to one file by master"
+cd write (iout,*) " before exchange at NSTEX intervals"
+cd endif
+ write (iout,*) "NREP= ",nrep
+ write (iout,*) "NSTEX= ",nstex
+ write (iout,*) "SYNC= ",mremdsync
+ write (iout,*) "NSYN= ",i_sync_step
+ write (iout,*) "TRAJCACHE= ",max_cache_traj_use
+ endif
+
+ t_exchange_only=(index(controlcard,'TONLY').gt.0)
+ call readi(controlcard,"HREMD",hremd,0)
+ if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
+ write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
+ endif
+ if(usampl.and.hremd.gt.0) then
+ write (iout,'(//a)')
+ & "========== ERROR: USAMPL and HREMD cannot be used together"
+#ifdef MPI
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#endif
+ stop
+ endif
+
+
+ remd_tlist=.false.
+ if (index(controlcard,'TLIST').gt.0) then
+ remd_tlist=.true.
+ call card_concat(controlcard1)
+ read(controlcard1,*) (remd_t(i),i=1,nrep)
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
+ endif
+ remd_mlist=.false.
+ if (index(controlcard,'MLIST').gt.0) then
+ remd_mlist=.true.
+ call card_concat(controlcard1)
+ read(controlcard1,*) (remd_m(i),i=1,nrep)
+ if(me.eq.king.or..not.out1file) then
+ write (iout,*)'mlist',(remd_m(i),i=1,nrep)
+ iremd_m_total=0
+ do i=1,nrep
+ iremd_m_total=iremd_m_total+remd_m(i)
+ enddo
+ if(hremd.gt.1)then
+ write (iout,*) 'Total number of replicas ',
+ & iremd_m_total*hremd
+ else
+ write (iout,*) 'Total number of replicas ',iremd_m_total
+ endif
+ endif
+ endif
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine read_MDpar
+C
+C Read MD settings
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.TIME1'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
+ include 'COMMON.GEO'
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ character*80 ucase
+ character*320 controlcard
+
+ call card_concat(controlcard)
+ call readi(controlcard,"NSTEP",n_timestep,1000000)
+ call readi(controlcard,"NTWE",ntwe,100)
+ call readi(controlcard,"NTWX",ntwx,1000)
+ call reada(controlcard,"DT",d_time,1.0d-1)
+ call reada(controlcard,"DVMAX",dvmax,2.0d1)
+ call reada(controlcard,"DAMAX",damax,1.0d1)
+ call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
+ call readi(controlcard,"LANG",lang,0)
+ RESPA = index(controlcard,"RESPA") .gt. 0
+ call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
+ ntime_split0=ntime_split
+ call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
+ ntime_split0=ntime_split
+ call reada(controlcard,"R_CUT",r_cut,2.0d0)
+ call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+ rest = index(controlcard,"REST").gt.0
+ tbf = index(controlcard,"TBF").gt.0
+ call readi(controlcard,"HMC",hmc,0)
+ tnp = index(controlcard,"NOSEPOINCARE99").gt.0
+ tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
+ tnh = index(controlcard,"NOSEHOOVER96").gt.0
+ if (RESPA.and.tnh)then
+ xiresp = index(controlcard,"XIRESP").gt.0
+ endif
+ call reada(controlcard,"Q_NP",Q_np,0.1d0)
+ usampl = index(controlcard,"USAMPL").gt.0
+
+ mdpdb = index(controlcard,"MDPDB").gt.0
+ call reada(controlcard,"T_BATH",t_bath,300.0d0)
+ call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
+ call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
+ call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
+ if (count_reset_moment.eq.0) count_reset_moment=1000000000
+ call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
+ reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
+ reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
+ if (count_reset_vel.eq.0) count_reset_vel=1000000000
+ large = index(controlcard,"LARGE").gt.0
+ print_compon = index(controlcard,"PRINT_COMPON").gt.0
+ rattle = index(controlcard,"RATTLE").gt.0
+c if performing umbrella sampling, fragments constrained are read from the fragment file
+ nset=0
+ if(usampl) then
+ call read_fragments
+ endif
+
+ if(me.eq.king.or..not.out1file) then
+ write (iout,*)
+ write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
+ write (iout,*)
+ write (iout,'(a)') "The units are:"
+ write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
+ write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
+ & " acceleration: angstrom/(48.9 fs)**2"
+ write (iout,'(a)') "energy: kcal/mol, temperature: K"
+ write (iout,*)
+ write (iout,'(a60,i10)') "Number of time steps:",n_timestep
+ write (iout,'(a60,f10.5,a)')
+ & "Initial time step of numerical integration:",d_time,
+ & " natural units"
+ write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
+ if (RESPA) then
+ write (iout,'(2a,i4,a)')
+ & "A-MTS algorithm used; initial time step for fast-varying",
+ & " short-range forces split into",ntime_split," steps."
+ write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
+ & r_cut," lambda",rlamb
+ endif
+ write (iout,'(2a,f10.5)')
+ & "Maximum acceleration threshold to reduce the time step",
+ & "/increase split number:",damax
+ write (iout,'(2a,f10.5)')
+ & "Maximum predicted energy drift to reduce the timestep",
+ & "/increase split number:",edriftmax
+ write (iout,'(a60,f10.5)')
+ & "Maximum velocity threshold to reduce velocities:",dvmax
+ write (iout,'(a60,i10)') "Frequency of property output:",ntwe
+ write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
+ if (rattle) write (iout,'(a60)')
+ & "Rattle algorithm used to constrain the virtual bonds"
+ endif
+ reset_fricmat=1000
+ if (lang.gt.0) then
+ call reada(controlcard,"ETAWAT",etawat,0.8904d0)
+ call reada(controlcard,"RWAT",rwat,1.4d0)
+ call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
+ surfarea=index(controlcard,"SURFAREA").gt.0
+ call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(/a,$)') "Langevin dynamics calculation"
+ if (lang.eq.1) then
+ write (iout,'(a/)')
+ & " with direct integration of Langevin equations"
+ else if (lang.eq.2) then
+ write (iout,'(a/)') " with TINKER stochasic MD integrator"
+ else if (lang.eq.3) then
+ write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
+ else if (lang.eq.4) then
+ write (iout,'(a/)') " in overdamped mode"
+ else
+ write (iout,'(//a,i5)')
+ & "=========== ERROR: Unknown Langevin dynamics mode:",lang
+ stop
+ endif
+ write (iout,'(a60,f10.5)') "Temperature:",t_bath
+ write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
+ write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
+ write (iout,'(a60,f10.5)')
+ & "Scaling factor of the friction forces:",scal_fric
+ if (surfarea) write (iout,'(2a,i10,a)')
+ & "Friction coefficients will be scaled by solvent-accessible",
+ & " surface area every",reset_fricmat," steps."
+ endif
+c Calculate friction coefficients and bounds of stochastic forces
+ eta=6*pi*cPoise*etawat
+ if(me.eq.king.or..not.out1file)
+ & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
+ & ,eta
+ gamp=scal_fric*(pstok+rwat)*eta
+ stdfp=dsqrt(2*Rb*t_bath/d_time)
+ do i=1,ntyp
+ gamsc(i)=scal_fric*(restok(i)+rwat)*eta
+ stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+ enddo
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(/2a/)')
+ & "Radii of site types and friction coefficients and std's of",
+ & " stochastic forces of fully exposed sites"
+ write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
+ do i=1,ntyp
+ write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
+ & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
+ enddo
+ endif
+ else if (tbf) then
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(a)') "Berendsen bath calculation"
+ write (iout,'(a60,f10.5)') "Temperature:",t_bath
+ write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
+ if (reset_moment)
+ & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
+ & count_reset_moment," steps"
+ if (reset_vel)
+ & write (iout,'(a,i10,a)')
+ & "Velocities will be reset at random every",count_reset_vel,
+ & " steps"
+ endif
+ else if (tnp .or. tnp1 .or. tnh) then
+ if (tnp .or. tnp1) then
+ write (iout,'(a)') "Nose-Poincare bath calculation"
+ if (tnp) write (iout,'(a)')
+ & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
+ if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
+ else
+ write (iout,'(a)') "Nose-Hoover bath calculation"
+ write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
+ nresn=1
+ nyosh=1
+ nnos=1
+ do i=1,nnos
+ qmass(i)=Q_np
+ xlogs(i)=1.0
+ vlogs(i)=0.0
+ enddo
+ do i=1,nyosh
+ WDTI(i) = 1.0*d_time/nresn
+ WDTI2(i)=WDTI(i)/2
+ WDTI4(i)=WDTI(i)/4
+ WDTI8(i)=WDTI(i)/8
+ enddo
+ if (RESPA) then
+ if(xiresp) then
+ write (iout,'(a)') "NVT-XI-RESPA algorithm"
+ else
+ write (iout,'(a)') "NVT-XO-RESPA algorithm"
+ endif
+ do i=1,nyosh
+ WDTIi(i) = 1.0*d_time/nresn/ntime_split
+ WDTIi2(i)=WDTIi(i)/2
+ WDTIi4(i)=WDTIi(i)/4
+ WDTIi8(i)=WDTIi(i)/8
+ enddo
+ endif
+ endif
+
+ write (iout,'(a60,f10.5)') "Temperature:",t_bath
+ write (iout,'(a60,f10.5)') "Q =",Q_np
+ if (reset_moment)
+ & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
+ & count_reset_moment," steps"
+ if (reset_vel)
+ & write (iout,'(a,i10,a)')
+ & "Velocities will be reset at random every",count_reset_vel,
+ & " steps"
+
+ else if (hmc.gt.0) then
+ write (iout,'(a)') "Hybrid Monte Carlo calculation"
+ write (iout,'(a60,f10.5)') "Temperature:",t_bath
+ write (iout,'(a60,i10)')
+ & "Number of MD steps between Metropolis tests:",hmc
+
+ else
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(a31)') "Microcanonical mode calculation"
+ endif
+ if(me.eq.king.or..not.out1file)then
+ if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
+ if (usampl) then
+ write(iout,*) "MD running with constraints."
+ write(iout,*) "Equilibration time ", eq_time, " mtus."
+ write(iout,*) "Constraining ", nfrag," fragments."
+ write(iout,*) "Length of each fragment, weight and q0:"
+ do iset=1,nset
+ write (iout,*) "Set of restraints #",iset
+ do i=1,nfrag
+ write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
+ & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
+ enddo
+ write(iout,*) "constraints between ", npair, "fragments."
+ write(iout,*) "constraint pairs, weights and q0:"
+ do i=1,npair
+ write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
+ & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
+ enddo
+ write(iout,*) "angle constraints within ", nfrag_back,
+ & "backbone fragments."
+ write(iout,*) "fragment, weights:"
+ do i=1,nfrag_back
+ write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
+ & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
+ & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
+ enddo
+ enddo
+ iset=mod(kolor,nset)+1
+ endif
+ endif
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine molread
+C
+C Read molecular data.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+ integer error_msg
+#endif
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.LOCAL'
+ include 'COMMON.NAMES'
+ include 'COMMON.CHAIN'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.CONTROL'
+ include 'COMMON.DBASE'
+ include 'COMMON.THREAD'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.TORCNSTR'
+ include 'COMMON.TIME1'
+ include 'COMMON.BOUNDS'
+ include 'COMMON.MD'
+ include 'COMMON.REMD'
+ include 'COMMON.SETUP'
+ character*4 sequence(maxres)
+ integer rescode
+ double precision x(maxvar)
+ character*256 pdbfile
+ character*320 weightcard
+ character*80 weightcard_t,ucase
+ dimension itype_pdb(maxres)
+ common /pizda/ itype_pdb
+ logical seq_comp,fail
+ double precision energia(0:n_ene)
+ integer ilen
+ external ilen
+C
+C Body
+C
+C Read weights of the subsequent energy terms.
+ if(hremd.gt.0) then
+
+ k=0
+ do il=1,hremd
+ do i=1,nrep
+ do j=1,remd_m(i)
+ i2set(k)=il
+ k=k+1
+ enddo
+ enddo
+ enddo
+
+ if(me.eq.king.or..not.out1file) then
+ write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
+ write (iout,*) 'Current weights for processor ',
+ & me,' set ',i2set(me)
+ endif
+
+ do i=1,hremd
+ call card_concat(weightcard)
+ call reada(weightcard,'WLONG',wlong,1.0D0)
+ call reada(weightcard,'WSC',wsc,wlong)
+ call reada(weightcard,'WSCP',wscp,wlong)
+ call reada(weightcard,'WELEC',welec,1.0D0)
+ call reada(weightcard,'WVDWPP',wvdwpp,welec)
+ call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
+ call reada(weightcard,'WCORR4',wcorr4,0.0D0)
+ call reada(weightcard,'WCORR5',wcorr5,0.0D0)
+ call reada(weightcard,'WCORR6',wcorr6,0.0D0)
+ call reada(weightcard,'WTURN3',wturn3,1.0D0)
+ call reada(weightcard,'WTURN4',wturn4,1.0D0)
+ call reada(weightcard,'WTURN6',wturn6,1.0D0)
+ call reada(weightcard,'WSCCOR',wsccor,1.0D0)
+ call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
+ call reada(weightcard,'WBOND',wbond,1.0D0)
+ call reada(weightcard,'WTOR',wtor,1.0D0)
+ call reada(weightcard,'WTORD',wtor_d,1.0D0)
+ call reada(weightcard,'WANG',wang,1.0D0)
+ call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+ call reada(weightcard,'SCAL14',scal14,0.4D0)
+ call reada(weightcard,'SCALSCP',scalscp,1.0d0)
+ call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
+ call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
+ call reada(weightcard,'TEMP0',temp0,300.0d0)
+ if (index(weightcard,'SOFT').gt.0) ipot=6
+C 12/1/95 Added weight for the multi-body term WCORR
+ call reada(weightcard,'WCORRH',wcorr,1.0D0)
+ if (wcorr4.gt.0.0d0) wcorr=wcorr4
+
+ hweights(i,1)=wsc
+ hweights(i,2)=wscp
+ hweights(i,3)=welec
+ hweights(i,4)=wcorr
+ hweights(i,5)=wcorr5
+ hweights(i,6)=wcorr6
+ hweights(i,7)=wel_loc
+ hweights(i,8)=wturn3
+ hweights(i,9)=wturn4
+ hweights(i,10)=wturn6
+ hweights(i,11)=wang
+ hweights(i,12)=wscloc
+ hweights(i,13)=wtor
+ hweights(i,14)=wtor_d
+ hweights(i,15)=wstrain
+ hweights(i,16)=wvdwpp
+ hweights(i,17)=wbond
+ hweights(i,18)=scal14
+ hweights(i,21)=wsccor
+
+ enddo
+
+ do i=1,n_ene
+ weights(i)=hweights(i2set(me),i)
+ enddo
+ wsc =weights(1)
+ wscp =weights(2)
+ welec =weights(3)
+ wcorr =weights(4)
+ wcorr5 =weights(5)
+ wcorr6 =weights(6)
+ wel_loc=weights(7)
+ wturn3 =weights(8)
+ wturn4 =weights(9)
+ wturn6 =weights(10)
+ wang =weights(11)
+ wscloc =weights(12)
+ wtor =weights(13)
+ wtor_d =weights(14)
+ wstrain=weights(15)
+ wvdwpp =weights(16)
+ wbond =weights(17)
+ scal14 =weights(18)
+ wsccor =weights(21)
+
+
+ else
+ call card_concat(weightcard)
+ call reada(weightcard,'WLONG',wlong,1.0D0)
+ call reada(weightcard,'WSC',wsc,wlong)
+ call reada(weightcard,'WSCP',wscp,wlong)
+ call reada(weightcard,'WELEC',welec,1.0D0)
+ call reada(weightcard,'WVDWPP',wvdwpp,welec)
+ call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
+ call reada(weightcard,'WCORR4',wcorr4,0.0D0)
+ call reada(weightcard,'WCORR5',wcorr5,0.0D0)
+ call reada(weightcard,'WCORR6',wcorr6,0.0D0)
+ call reada(weightcard,'WTURN3',wturn3,1.0D0)
+ call reada(weightcard,'WTURN4',wturn4,1.0D0)
+ call reada(weightcard,'WTURN6',wturn6,1.0D0)
+ call reada(weightcard,'WSCCOR',wsccor,1.0D0)
+ call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
+ call reada(weightcard,'WBOND',wbond,1.0D0)
+ call reada(weightcard,'WTOR',wtor,1.0D0)
+ call reada(weightcard,'WTORD',wtor_d,1.0D0)
+ call reada(weightcard,'WANG',wang,1.0D0)
+ call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+ call reada(weightcard,'SCAL14',scal14,0.4D0)
+ call reada(weightcard,'SCALSCP',scalscp,1.0d0)
+ call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
+ call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
+ call reada(weightcard,'TEMP0',temp0,300.0d0)
+ if (index(weightcard,'SOFT').gt.0) ipot=6
+C 12/1/95 Added weight for the multi-body term WCORR
+ call reada(weightcard,'WCORRH',wcorr,1.0D0)
+ if (wcorr4.gt.0.0d0) wcorr=wcorr4
+ weights(1)=wsc
+ weights(2)=wscp
+ weights(3)=welec
+ weights(4)=wcorr
+ weights(5)=wcorr5
+ weights(6)=wcorr6
+ weights(7)=wel_loc
+ weights(8)=wturn3
+ weights(9)=wturn4
+ weights(10)=wturn6
+ weights(11)=wang
+ weights(12)=wscloc
+ weights(13)=wtor
+ weights(14)=wtor_d
+ weights(15)=wstrain
+ weights(16)=wvdwpp
+ weights(17)=wbond
+ weights(18)=scal14
+ weights(21)=wsccor
+ endif
+
+ if(me.eq.king.or..not.out1file)
+ & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
+ & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
+ & wturn4,wturn6
+ 10 format (/'Energy-term weights (unscaled):'//
+ & 'WSCC= ',f10.6,' (SC-SC)'/
+ & 'WSCP= ',f10.6,' (SC-p)'/
+ & 'WELEC= ',f10.6,' (p-p electr)'/
+ & 'WVDWPP= ',f10.6,' (p-p VDW)'/
+ & 'WBOND= ',f10.6,' (stretching)'/
+ & 'WANG= ',f10.6,' (bending)'/
+ & 'WSCLOC= ',f10.6,' (SC local)'/
+ & 'WTOR= ',f10.6,' (torsional)'/
+ & 'WTORD= ',f10.6,' (double torsional)'/
+ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
+ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
+ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
+ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
+ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
+ & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
+ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
+ & 'WTURN4= ',f10.6,' (turns, 4th order)'/
+ & 'WTURN6= ',f10.6,' (turns, 6th order)')
+ if(me.eq.king.or..not.out1file)then
+ if (wcorr4.gt.0.0d0) then
+ write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
+ & 'between contact pairs of peptide groups'
+ write (iout,'(2(a,f5.3/))')
+ & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
+ & 'Range of quenching the correlation terms:',2*delt_corr
+ else if (wcorr.gt.0.0d0) then
+ write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
+ & 'between contact pairs of peptide groups'
+ endif
+ write (iout,'(a,f8.3)')
+ & 'Scaling factor of 1,4 SC-p interactions:',scal14
+ write (iout,'(a,f8.3)')
+ & 'General scaling factor of SC-p interactions:',scalscp
+ endif
+ r0_corr=cutoff_corr-delt_corr
+ do i=1,20
+ aad(i,1)=scalscp*aad(i,1)
+ aad(i,2)=scalscp*aad(i,2)
+ bad(i,1)=scalscp*bad(i,1)
+ bad(i,2)=scalscp*bad(i,2)
+ enddo
+ call rescale_weights(t_bath)
+ if(me.eq.king.or..not.out1file)
+ & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
+ & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
+ & wturn4,wturn6
+ 22 format (/'Energy-term weights (scaled):'//
+ & 'WSCC= ',f10.6,' (SC-SC)'/
+ & 'WSCP= ',f10.6,' (SC-p)'/
+ & 'WELEC= ',f10.6,' (p-p electr)'/
+ & 'WVDWPP= ',f10.6,' (p-p VDW)'/
+ & 'WBOND= ',f10.6,' (stretching)'/
+ & 'WANG= ',f10.6,' (bending)'/
+ & 'WSCLOC= ',f10.6,' (SC local)'/
+ & 'WTOR= ',f10.6,' (torsional)'/
+ & 'WTORD= ',f10.6,' (double torsional)'/
+ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
+ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
+ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
+ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
+ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
+ & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
+ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
+ & 'WTURN4= ',f10.6,' (turns, 4th order)'/
+ & 'WTURN6= ',f10.6,' (turns, 6th order)')
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) "Reference temperature for weights calculation:",
+ & temp0
+ call reada(weightcard,"D0CM",d0cm,3.78d0)
+ call reada(weightcard,"AKCM",akcm,15.1d0)
+ call reada(weightcard,"AKTH",akth,11.0d0)
+ call reada(weightcard,"AKCT",akct,12.0d0)
+ call reada(weightcard,"V1SS",v1ss,-1.08d0)
+ call reada(weightcard,"V2SS",v2ss,7.61d0)
+ call reada(weightcard,"V3SS",v3ss,13.7d0)
+ call reada(weightcard,"EBR",ebr,-5.50D0)
+ if(me.eq.king.or..not.out1file) then
+ write (iout,*) "Parameters of the SS-bond potential:"
+ write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
+ & " AKCT",akct
+ write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
+ write (iout,*) "EBR",ebr
+ print *,'indpdb=',indpdb,' pdbref=',pdbref
+ endif
+ if (indpdb.gt.0 .or. pdbref) then
+ read(inp,'(a)') pdbfile
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(2a)') 'PDB data will be read from file ',
+ & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a)') 'Error opening PDB file.'
+ stop
+ 34 continue
+c print *,'Begin reading pdb data'
+ call readpdb
+c print *,'Finished reading pdb data'
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
+ & ' nstart_sup=',nstart_sup
+ do i=1,nres
+ itype_pdb(i)=itype(i)
+ enddo
+ close (ipdbin)
+ nnt=nstart_sup
+ nct=nstart_sup+nsup-1
+ call contact(.false.,ncont_ref,icont_ref,co)
+
+ if (sideadd) then
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*)'Adding sidechains'
+ maxsi=1000
+ do i=2,nres-1
+ iti=itype(i)
+ if (iti.ne.10) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+ call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
+ nsi=nsi+1
+ enddo
+ if(fail) write(iout,*)'Adding sidechain failed for res ',
+ & i,' after ',nsi,' trials'
+ endif
+ enddo
+ endif
+ endif
+ if (indpdb.eq.0) then
+C Read sequence if not taken from the pdb file.
+ read (inp,*) nres
+c print *,'nres=',nres
+ if (iscode.gt.0) then
+ read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
+ else
+ read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
+ endif
+C Convert sequence to numeric code
+ do i=1,nres
+ itype(i)=rescode(i,sequence(i),iscode)
+ enddo
+C Assign initial virtual bond lengths
+ do i=2,nres
+ vbld(i)=vbl
+ vbld_inv(i)=vblinv
+ enddo
+ do i=2,nres-1
+ vbld(i+nres)=dsc(itype(i))
+ vbld_inv(i+nres)=dsc_inv(itype(i))
+c write (iout,*) "i",i," itype",itype(i),
+c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
+ enddo
+ endif
+c print *,nres
+c print '(20i4)',(itype(i),i=1,nres)
+ do i=1,nres
+#ifdef PROCOR
+ if (itype(i).eq.21 .or. itype(i+1).eq.21) then
+#else
+ if (itype(i).eq.21) then
+#endif
+ itel(i)=0
+#ifdef PROCOR
+ else if (itype(i+1).ne.20) then
+#else
+ else if (itype(i).ne.20) then
+#endif
+ itel(i)=1
+ else
+ itel(i)=2
+ endif
+ enddo
+ if(me.eq.king.or..not.out1file)then
+ write (iout,*) "ITEL"
+ do i=1,nres-1
+ write (iout,*) i,itype(i),itel(i)
+ enddo
+ print *,'Call Read_Bridge.'
+ endif
+ call read_bridge
+C 8/13/98 Set limits to generating the dihedral angles
+ do i=1,nres
+ phibound(1,i)=-pi
+ phibound(2,i)=pi
+ enddo
+ read (inp,*) ndih_constr
+ if (ndih_constr.gt.0) then
+ read (inp,*) ftors
+ read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
+ if(me.eq.king.or..not.out1file)then
+ write (iout,*)
+ & 'There are',ndih_constr,' constraints on phi angles.'
+ do i=1,ndih_constr
+ write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
+ enddo
+ endif
+ do i=1,ndih_constr
+ phi0(i)=deg2rad*phi0(i)
+ drange(i)=deg2rad*drange(i)
+ enddo
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) 'FTORS',ftors
+ do i=1,ndih_constr
+ ii = idih_constr(i)
+ phibound(1,ii) = phi0(i)-drange(i)
+ phibound(2,ii) = phi0(i)+drange(i)
+ enddo
+ endif
+ nnt=1
+#ifdef MPI
+ if (me.eq.king) then
+#endif
+ write (iout,'(a)') 'Boundaries in phi angle sampling:'
+ do i=1,nres
+ write (iout,'(a3,i5,2f10.1)')
+ & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
+ enddo
+#ifdef MP
+ endif
+#endif
+ nct=nres
+cd print *,'NNT=',NNT,' NCT=',NCT
+ if (itype(1).eq.21) nnt=2
+ if (itype(nres).eq.21) nct=nct-1
+ if (pdbref) then
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(a,i3)') 'nsup=',nsup
+ nstart_seq=nnt
+ if (nsup.le.(nct-nnt+1)) then
+ do i=0,nct-nnt+1-nsup
+ if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
+ nstart_seq=nnt+i
+ goto 111
+ endif
+ enddo
+ write (iout,'(a)')
+ & 'Error - sequences to be superposed do not match.'
+ stop
+ else
+ do i=0,nsup-(nct-nnt+1)
+ if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
+ & then
+ nstart_sup=nstart_sup+i
+ nsup=nct-nnt+1
+ goto 111
+ endif
+ enddo
+ write (iout,'(a)')
+ & 'Error - sequences to be superposed do not match.'
+ endif
+ 111 continue
+ if (nsup.eq.0) nsup=nct-nnt
+ if (nstart_sup.eq.0) nstart_sup=nnt
+ if (nstart_seq.eq.0) nstart_seq=nnt
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
+ & ' nstart_seq=',nstart_seq
+ endif
+c--- Zscore rms -------
+ if (nz_start.eq.0) nz_start=nnt
+ if (nz_end.eq.0 .and. nsup.gt.0) then
+ nz_end=nnt+nsup-1
+ else if (nz_end.eq.0) then
+ nz_end=nct
+ endif
+ if(me.eq.king.or..not.out1file)then
+ write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
+ write (iout,*) 'IZ_SC=',iz_sc
+ endif
+c----------------------
+ call init_int_table
+ if (refstr) then
+ if (.not.pdbref) then
+ call read_angles(inp,*38)
+ goto 39
+ 38 write (iout,'(a)') 'Error reading reference structure.'
+#ifdef MPI
+ call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+ stop 'Error reading reference structure'
+#endif
+ 39 call chainbuild
+ call setup_var
+czscore call geom_to_var(nvar,coord_exp_zs(1,1))
+ nstart_sup=nnt
+ nstart_seq=nnt
+ nsup=nct-nnt+1
+ do i=1,2*nres
+ do j=1,3
+ cref(j,i)=c(j,i)
+ enddo
+ enddo
+ call contact(.true.,ncont_ref,icont_ref,co)
+ endif
+c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+ call flush(iout)
+ if (constr_dist.gt.0) call read_dist_constr
+c write (iout,*) "After read_dist_constr nhpb",nhpb
+ call hpb_partition
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) 'Contact order:',co
+ if (pdbref) then
+ if(me.eq.king.or..not.out1file)
+ & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
+ do i=1,ncont_ref
+ do j=1,2
+ icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
+ enddo
+ if(me.eq.king.or..not.out1file)
+ & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
+ & icont_ref(1,i),' ',
+ & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
+ enddo
+ endif
+ endif
+ if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
+ & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
+ & modecalc.ne.10) then
+C If input structure hasn't been supplied from the PDB file read or generate
+C initial geometry.
+ if (iranconf.eq.0 .and. .not. extconf) then
+ if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
+ & write (iout,'(a)') 'Initial geometry will be read in.'
+ if (read_cart) then
+ read(inp,'(8f10.5)',end=36,err=36)
+ & ((c(l,k),l=1,3),k=1,nres),
+ & ((c(l,k+nres),l=1,3),k=nnt,nct)
+ call int_from_cart1(.false.)
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ endif
+ enddo
+ return
+ else
+ call read_angles(inp,*36)
+ endif
+ goto 37
+ 36 write (iout,'(a)') 'Error reading angle file.'
+#ifdef MPI
+ call mpi_finalize( MPI_COMM_WORLD,IERR )
+#endif
+ stop 'Error reading angle file.'
+ 37 continue
+ else if (extconf) then
+ if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
+ & write (iout,'(a)') 'Extended chain initial geometry.'
+ do i=3,nres
+ theta(i)=90d0*deg2rad
+ enddo
+ do i=4,nres
+ phi(i)=180d0*deg2rad
+ enddo
+ do i=2,nres-1
+ alph(i)=110d0*deg2rad
+ enddo
+ do i=2,nres-1
+ omeg(i)=-120d0*deg2rad
+ enddo
+ else
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(a)') 'Random-generated initial geometry.'
+
+
+#ifdef MPI
+ if (me.eq.king .or. fg_rank.eq.0 .and. (
+ & modecalc.eq.12 .or. modecalc.eq.14) ) then
+#endif
+ do itrial=1,100
+ itmp=1
+ call gen_rand_conf(itmp,*30)
+ goto 40
+ 30 write (iout,*) 'Failed to generate random conformation',
+ & ', itrial=',itrial
+ write (*,*) 'Processor:',me,
+ & ' Failed to generate random conformation',
+ & ' itrial=',itrial
+ call intout
+
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+ enddo
+ write (iout,'(a,i3,a)') 'Processor:',me,
+ & ' error in generating random conformation.'
+ write (*,'(a,i3,a)') 'Processor:',me,
+ & ' error in generating random conformation.'
+ call flush(iout)
+#ifdef MPI
+ call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+ 40 continue
+ endif
+#else
+ 40 continue
+#endif
+ endif
+ elseif (modecalc.eq.4) then
+ read (inp,'(a)') intinname
+ open (intin,file=intinname,status='old',err=333)
+ if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
+ & write (iout,'(a)') 'intinname',intinname
+ write (*,'(a)') 'Processor',myrank,' intinname',intinname
+ goto 334
+ 333 write (iout,'(2a)') 'Error opening angle file ',intinname
+#ifdef MPI
+ call MPI_Finalize(MPI_COMM_WORLD,IERR)
+#endif
+ stop 'Error opening angle file.'
+ 334 continue
+
+ endif
+C Generate distance constraints, if the PDB structure is to be regularized.
+ if (nthread.gt.0) then
+ call read_threadbase
+ endif
+ call setup_var
+ if (me.eq.king .or. .not. out1file)
+ & call intout
+ if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
+ write (iout,'(/a,i3,a)')
+ & 'The chain contains',ns,' disulfide-bridging cysteines.'
+ write (iout,'(20i4)') (iss(i),i=1,ns)
+ write (iout,'(/a/)') 'Pre-formed links are:'
+ do i=1,nss
+ i1=ihpb(i)-nres
+ i2=jhpb(i)-nres
+ it1=itype(i1)
+ it2=itype(i2)
+ if (me.eq.king.or..not.out1file)
+ & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
+ & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
+ & ebr,forcon(i)
+ enddo
+ write (iout,'(a)')
+ endif
+ if (i2ndstr.gt.0) call secstrp2dihc
+c call geom_to_var(nvar,x)
+c call etotal(energia(0))
+c call enerprint(energia(0))
+c call briefout(0,etot)
+c stop
+cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
+cd write (iout,'(a)') 'Variable list:'
+cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
+#ifdef MPI
+ if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
+ & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
+ & 'Processor',myrank,': end reading molecular data.'
+#endif
+ return
+ end
+c--------------------------------------------------------------------------
+ logical function seq_comp(itypea,itypeb,length)
+ implicit none
+ integer length,itypea(length),itypeb(length)
+ integer i
+ do i=1,length
+ if (itypea(i).ne.itypeb(i)) then
+ seq_comp=.false.
+ return
+ endif
+ enddo
+ seq_comp=.true.
+ return
+ end
+c-----------------------------------------------------------------------------
+ subroutine read_bridge
+C Read information about disulfide bridges.
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.LOCAL'
+ include 'COMMON.NAMES'
+ include 'COMMON.CHAIN'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.CONTROL'
+ include 'COMMON.DBASE'
+ include 'COMMON.THREAD'
+ include 'COMMON.TIME1'
+ include 'COMMON.SETUP'
+C Read bridging residues.
+ read (inp,*) ns,(iss(i),i=1,ns)
+ print *,'ns=',ns
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
+C Check whether the specified bridging residues are cystines.
+ do i=1,ns
+ if (itype(iss(i)).ne.1) then
+ if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
+ & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
+ & ' can form a disulfide bridge?!!!'
+ write (*,'(2a,i3,a)')
+ & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
+ & ' can form a disulfide bridge?!!!'
+#ifdef MPI
+ call MPI_Finalize(MPI_COMM_WORLD,ierror)
+ stop
+#endif
+ endif
+ enddo
+C Read preformed bridges.
+ if (ns.gt.0) then
+ read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
+ write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
+ if (nss.gt.0) then
+ nhpb=nss
+C Check if the residues involved in bridges are in the specified list of
+C bridging residues.
+ do i=1,nss
+ do j=1,i-1
+ if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
+ & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
+ write (iout,'(a,i3,a)') 'Disulfide pair',i,
+ & ' contains residues present in other pairs.'
+ write (*,'(a,i3,a)') 'Disulfide pair',i,
+ & ' contains residues present in other pairs.'
+#ifdef MPI
+ call MPI_Finalize(MPI_COMM_WORLD,ierror)
+ stop
+#endif
+ endif
+ enddo
+ do j=1,ns
+ if (ihpb(i).eq.iss(j)) goto 10
+ enddo
+ write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
+ 10 continue
+ do j=1,ns
+ if (jhpb(i).eq.iss(j)) goto 20
+ enddo
+ write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
+ 20 continue
+ dhpb(i)=dbr
+ forcon(i)=fbr
+ enddo
+ do i=1,nss
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ enddo
+ endif
+ endif
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine read_x(kanal,*)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+c Read coordinates from input
+c
+ read(kanal,'(8f10.5)',end=10,err=10)
+ & ((c(l,k),l=1,3),k=1,nres),
+ & ((c(l,k+nres),l=1,3),k=nnt,nct)
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ enddo
+ call int_from_cart1(.false.)
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ endif
+ enddo
+
+ return
+ 10 return1
+ end
+c----------------------------------------------------------------------------
+ subroutine read_threadbase
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.LOCAL'
+ include 'COMMON.NAMES'
+ include 'COMMON.CHAIN'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.CONTROL'
+ include 'COMMON.DBASE'
+ include 'COMMON.THREAD'
+ include 'COMMON.TIME1'
+C Read pattern database for threading.
+ read (icbase,*) nseq
+ do i=1,nseq
+ read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
+ & nres_base(2,i),nres_base(3,i)
+ read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
+ & nres_base(1,i))
+c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
+c & nres_base(2,i),nres_base(3,i)
+c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
+c & nres_base(1,i))
+ enddo
+ close (icbase)
+ if (weidis.eq.0.0D0) weidis=0.1D0
+ do i=nnt,nct
+ do j=i+2,nct
+ nhpb=nhpb+1
+ ihpb(nhpb)=i
+ jhpb(nhpb)=j
+ forcon(nhpb)=weidis
+ enddo
+ enddo
+ read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
+ write (iout,'(a,i5)') 'nexcl: ',nexcl
+ write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine setup_var
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.LOCAL'
+ include 'COMMON.NAMES'
+ include 'COMMON.CHAIN'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.CONTROL'
+ include 'COMMON.DBASE'
+ include 'COMMON.THREAD'
+ include 'COMMON.TIME1'
+C Set up variable list.
+ ntheta=nres-2
+ nphi=nres-3
+ nvar=ntheta+nphi
+ nside=0
+ do i=2,nres-1
+ if (itype(i).ne.10) then
+ nside=nside+1
+ ialph(i,1)=nvar+nside
+ ialph(nside,2)=i
+ endif
+ enddo
+ if (indphi.gt.0) then
+ nvar=nphi
+ else if (indback.gt.0) then
+ nvar=nphi+ntheta
+ else
+ nvar=nvar+2*nside
+ endif
+cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine gen_dist_constr
+C Generate CA distance constraints.
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.LOCAL'
+ include 'COMMON.NAMES'
+ include 'COMMON.CHAIN'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.HEADER'
+ include 'COMMON.CONTROL'
+ include 'COMMON.DBASE'
+ include 'COMMON.THREAD'
+ include 'COMMON.TIME1'
+ dimension itype_pdb(maxres)
+ common /pizda/ itype_pdb
+ character*2 iden
+cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
+cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
+cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
+cd & ' nsup',nsup
+ do i=nstart_sup,nstart_sup+nsup-1
+cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
+cd & ' seq_pdb', restyp(itype_pdb(i))
+ do j=i+2,nstart_sup+nsup-1
+ nhpb=nhpb+1
+ ihpb(nhpb)=i+nstart_seq-nstart_sup
+ jhpb(nhpb)=j+nstart_seq-nstart_sup
+ forcon(nhpb)=weidis
+ dhpb(nhpb)=dist(i,j)
+ enddo
+ enddo
+cd write (iout,'(a)') 'Distance constraints:'
+cd do i=nss+1,nhpb
+cd ii=ihpb(i)
+cd jj=jhpb(i)
+cd iden='CA'
+cd if (ii.gt.nres) then
+cd iden='SC'
+cd ii=ii-nres
+cd jj=jj-nres
+cd endif
+cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
+cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
+cd & dhpb(i),forcon(i)
+cd enddo
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine map_read
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MAP'
+ include 'COMMON.IOUNITS'
+ character*3 angid(4) /'THE','PHI','ALP','OME'/
+ character*80 mapcard,ucase
+ do imap=1,nmap
+ read (inp,'(a)') mapcard
+ mapcard=ucase(mapcard)
+ if (index(mapcard,'PHI').gt.0) then
+ kang(imap)=1
+ else if (index(mapcard,'THE').gt.0) then
+ kang(imap)=2
+ else if (index(mapcard,'ALP').gt.0) then
+ kang(imap)=3
+ else if (index(mapcard,'OME').gt.0) then
+ kang(imap)=4
+ else
+ write(iout,'(a)')'Error - illegal variable spec in MAP card.'
+ stop 'Error - illegal variable spec in MAP card.'
+ endif
+ call readi (mapcard,'RES1',res1(imap),0)
+ call readi (mapcard,'RES2',res2(imap),0)
+ if (res1(imap).eq.0) then
+ res1(imap)=res2(imap)
+ else if (res2(imap).eq.0) then
+ res2(imap)=res1(imap)
+ endif
+ if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
+ write (iout,'(a)')
+ & 'Error - illegal definition of variable group in MAP.'
+ stop 'Error - illegal definition of variable group in MAP.'
+ endif
+ call reada(mapcard,'FROM',ang_from(imap),0.0D0)
+ call reada(mapcard,'TO',ang_to(imap),0.0D0)
+ call readi(mapcard,'NSTEP',nstep(imap),0)
+ if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
+ write (iout,'(a)')
+ & 'Illegal boundary and/or step size specification in MAP.'
+ stop 'Illegal boundary and/or step size specification in MAP.'
+ endif
+ enddo ! imap
+ return
+ end
+c----------------------------------------------------------------------------
+csa subroutine csaread
+csa implicit real*8 (a-h,o-z)
+csa include 'DIMENSIONS'
+csa include 'COMMON.IOUNITS'
+csa include 'COMMON.GEO'
+csa include 'COMMON.CSA'
+csa include 'COMMON.BANK'
+csa include 'COMMON.CONTROL'
+csa character*80 ucase
+csa character*620 mcmcard
+csa call card_concat(mcmcard)
+csa
+csa call readi(mcmcard,'NCONF',nconf,50)
+csa call readi(mcmcard,'NADD',nadd,0)
+csa call readi(mcmcard,'JSTART',jstart,1)
+csa call readi(mcmcard,'JEND',jend,1)
+csa call readi(mcmcard,'NSTMAX',nstmax,500000)
+csa call readi(mcmcard,'N0',n0,1)
+csa call readi(mcmcard,'N1',n1,6)
+csa call readi(mcmcard,'N2',n2,4)
+csa call readi(mcmcard,'N3',n3,0)
+csa call readi(mcmcard,'N4',n4,0)
+csa call readi(mcmcard,'N5',n5,0)
+csa call readi(mcmcard,'N6',n6,10)
+csa call readi(mcmcard,'N7',n7,0)
+csa call readi(mcmcard,'N8',n8,0)
+csa call readi(mcmcard,'N9',n9,0)
+csa call readi(mcmcard,'N14',n14,0)
+csa call readi(mcmcard,'N15',n15,0)
+csa call readi(mcmcard,'N16',n16,0)
+csa call readi(mcmcard,'N17',n17,0)
+csa call readi(mcmcard,'N18',n18,0)
+csa
+csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
+csa
+csa call readi(mcmcard,'NDIFF',ndiff,2)
+csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
+csa call readi(mcmcard,'IS1',is1,1)
+csa call readi(mcmcard,'IS2',is2,8)
+csa call readi(mcmcard,'NRAN0',nran0,4)
+csa call readi(mcmcard,'NRAN1',nran1,2)
+csa call readi(mcmcard,'IRR',irr,1)
+csa call readi(mcmcard,'NSEED',nseed,20)
+csa call readi(mcmcard,'NTOTAL',ntotal,10000)
+csa call reada(mcmcard,'CUT1',cut1,2.0d0)
+csa call reada(mcmcard,'CUT2',cut2,5.0d0)
+csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
+csa call readi(mcmcard,'ICMAX',icmax,3)
+csa call readi(mcmcard,'IRESTART',irestart,0)
+csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
+csa ntbankm=0
+csac!bankt
+csa call reada(mcmcard,'DELE',dele,20.0d0)
+csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
+csa call readi(mcmcard,'IREF',iref,0)
+csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
+csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
+csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
+csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
+csa write (iout,*) "NCONF_IN",nconf_in
+csa return
+csa end
+c----------------------------------------------------------------------------
+cfmc subroutine mcmfread
+cfmc implicit real*8 (a-h,o-z)
+cfmc include 'DIMENSIONS'
+cfmc include 'COMMON.MCMF'
+cfmc include 'COMMON.IOUNITS'
+cfmc include 'COMMON.GEO'
+cfmc character*80 ucase
+cfmc character*620 mcmcard
+cfmc call card_concat(mcmcard)
+cfmc
+cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
+cfmc write(iout,*)'MAXRANT=',maxrant
+cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
+cfmc write(iout,*)'MAXFAM=',maxfam
+cfmc call readi(mcmcard,'NNET1',nnet1,5)
+cfmc write(iout,*)'NNET1=',nnet1
+cfmc call readi(mcmcard,'NNET2',nnet2,4)
+cfmc write(iout,*)'NNET2=',nnet2
+cfmc call readi(mcmcard,'NNET3',nnet3,4)
+cfmc write(iout,*)'NNET3=',nnet3
+cfmc call readi(mcmcard,'ILASTT',ilastt,0)
+cfmc write(iout,*)'ILASTT=',ilastt
+cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
+cfmc write(iout,*)'MAXSTR=',maxstr
+cfmc maxstr_f=maxstr/maxfam
+cfmc write(iout,*)'MAXSTR_F=',maxstr_f
+cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
+cfmc write(iout,*)'NMCMF=',nmcmf
+cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
+cfmc write(iout,*)'IFOCUS=',ifocus
+cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
+cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
+cfmc call readi(mcmcard,'INTPRT',intprt,1000)
+cfmc write(iout,*)'INTPRT=',intprt
+cfmc call readi(mcmcard,'IPRT',iprt,100)
+cfmc write(iout,*)'IPRT=',iprt
+cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
+cfmc write(iout,*)'IMAXTR=',imaxtr
+cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
+cfmc write(iout,*)'MAXEVEN=',maxeven
+cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
+cfmc write(iout,*)'MAXEVEN1=',maxeven1
+cfmc call readi(mcmcard,'INIMIN',inimin,200)
+cfmc write(iout,*)'INIMIN=',inimin
+cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
+cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
+cfmc call readi(mcmcard,'NTHREAD',nthread,5)
+cfmc write(iout,*)'NTHREAD=',nthread
+cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
+cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
+cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
+cfmc write(iout,*)'MAXPERT=',maxpert
+cfmc call readi(mcmcard,'IRMSD',irmsd,1)
+cfmc write(iout,*)'IRMSD=',irmsd
+cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
+cfmc write(iout,*)'DENEMIN=',denemin
+cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
+cfmc write(iout,*)'RCUT1S=',rcut1s
+cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
+cfmc write(iout,*)'RCUT1E=',rcut1e
+cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
+cfmc write(iout,*)'RCUT2S=',rcut2s
+cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
+cfmc write(iout,*)'RCUT2E=',rcut2e
+cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
+cfmc write(iout,*)'DPERT1=',d_pert1
+cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
+cfmc write(iout,*)'DPERT1A=',d_pert1a
+cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
+cfmc write(iout,*)'DPERT2=',d_pert2
+cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
+cfmc write(iout,*)'DPERT2A=',d_pert2a
+cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
+cfmc write(iout,*)'DPERT2B=',d_pert2b
+cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
+cfmc write(iout,*)'DPERT2C=',d_pert2c
+cfmc d_pert1=deg2rad*d_pert1
+cfmc d_pert1a=deg2rad*d_pert1a
+cfmc d_pert2=deg2rad*d_pert2
+cfmc d_pert2a=deg2rad*d_pert2a
+cfmc d_pert2b=deg2rad*d_pert2b
+cfmc d_pert2c=deg2rad*d_pert2c
+cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
+cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
+cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
+cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
+cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
+cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
+cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
+cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
+cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
+cfmc write(iout,*)'RCUTINI=',rcutini
+cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
+cfmc write(iout,*)'GRAT=',grat
+cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
+cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
+cfmc
+cfmc return
+cfmc end
+c----------------------------------------------------------------------------
+ subroutine mcmread
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MCM'
+ include 'COMMON.MCE'
+ include 'COMMON.IOUNITS'
+ character*80 ucase
+ character*320 mcmcard
+ call card_concat(mcmcard)
+ call readi(mcmcard,'MAXACC',maxacc,100)
+ call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
+ call readi(mcmcard,'MAXTRIAL',maxtrial,100)
+ call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
+ call readi(mcmcard,'MAXREPM',maxrepm,200)
+ call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
+ call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
+ call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
+ call reada(mcmcard,'E_UP',e_up,5.0D0)
+ call reada(mcmcard,'DELTE',delte,0.1D0)
+ call readi(mcmcard,'NSWEEP',nsweep,5)
+ call readi(mcmcard,'NSTEPH',nsteph,0)
+ call readi(mcmcard,'NSTEPC',nstepc,0)
+ call reada(mcmcard,'TMIN',tmin,298.0D0)
+ call reada(mcmcard,'TMAX',tmax,298.0D0)
+ call readi(mcmcard,'NWINDOW',nwindow,0)
+ call readi(mcmcard,'PRINT_MC',print_mc,0)
+ print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
+ print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
+ ent_read=(index(mcmcard,'ENT_READ').gt.0)
+ call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
+ call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
+ call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
+ call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
+ call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
+ if (nwindow.gt.0) then
+ read (inp,*) (winstart(i),winend(i),i=1,nwindow)
+ do i=1,nwindow
+ winlen(i)=winend(i)-winstart(i)+1
+ enddo
+ endif
+ if (tmax.lt.tmin) tmax=tmin
+ if (tmax.eq.tmin) then
+ nstepc=0
+ nsteph=0
+ endif
+ if (nstepc.gt.0 .and. nsteph.gt.0) then
+ tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
+ tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
+ endif
+C Probabilities of different move types
+ sumpro_type(0)=0.0D0
+ call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
+ call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
+ sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
+ call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
+ sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
+ call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
+ sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
+ do i=1,MaxMoveType
+ print *,'i',i,' sumprotype',sumpro_type(i)
+ sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
+ print *,'i',i,' sumprotype',sumpro_type(i)
+ enddo
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine read_minim
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MINIM'
+ include 'COMMON.IOUNITS'
+ character*80 ucase
+ character*320 minimcard
+ call card_concat(minimcard)
+ call readi(minimcard,'MAXMIN',maxmin,2000)
+ call readi(minimcard,'MAXFUN',maxfun,5000)
+ call readi(minimcard,'MINMIN',minmin,maxmin)
+ call readi(minimcard,'MINFUN',minfun,maxmin)
+ call reada(minimcard,'TOLF',tolf,1.0D-2)
+ call reada(minimcard,'RTOLF',rtolf,1.0D-4)
+ write (iout,'(/80(1h*)/20x,a/80(1h*))')
+ & 'Options in energy minimization:'
+ write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
+ & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
+ & 'MinMin:',MinMin,' MinFun:',MinFun,
+ & ' TolF:',TolF,' RTolF:',RTolF
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine read_angles(kanal,*)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
+c Read angles from input
+c
+ read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
+ read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
+ read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
+ read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
+
+ do i=1,nres
+c 9/7/01 avoid 180 deg valence angle
+ if (theta(i).gt.179.99d0) theta(i)=179.99d0
+c
+ theta(i)=deg2rad*theta(i)
+ phi(i)=deg2rad*phi(i)
+ alph(i)=deg2rad*alph(i)
+ omeg(i)=deg2rad*omeg(i)
+ enddo
+ return
+ 10 return1
+ end
+c----------------------------------------------------------------------------
+ subroutine reada(rekord,lancuch,wartosc,default)
+ implicit none
+ character*(*) rekord,lancuch
+ double precision wartosc,default
+ integer ilen,iread
+ external ilen
+ iread=index(rekord,lancuch)
+ if (iread.eq.0) then
+ wartosc=default
+ return
+ endif
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,err=10,end=10) wartosc
+ return
+ 10 wartosc=default
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine readi(rekord,lancuch,wartosc,default)
+ implicit none
+ character*(*) rekord,lancuch
+ integer wartosc,default
+ integer ilen,iread
+ external ilen
+ iread=index(rekord,lancuch)
+ if (iread.eq.0) then
+ wartosc=default
+ return
+ endif
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,err=10,end=10) wartosc
+ return
+ 10 wartosc=default
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine multreadi(rekord,lancuch,tablica,dim,default)
+ implicit none
+ integer dim,i
+ integer tablica(dim),default
+ character*(*) rekord,lancuch
+ character*80 aux
+ integer ilen,iread
+ external ilen
+ do i=1,dim
+ tablica(i)=default
+ enddo
+ iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+ if (iread.eq.0) return
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+ 10 return
+ end
+c----------------------------------------------------------------------------
+ subroutine multreada(rekord,lancuch,tablica,dim,default)
+ implicit none
+ integer dim,i
+ double precision tablica(dim),default
+ character*(*) rekord,lancuch
+ character*80 aux
+ integer ilen,iread
+ external ilen
+ do i=1,dim
+ tablica(i)=default
+ enddo
+ iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+ if (iread.eq.0) return
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+ 10 return
+ end
+c----------------------------------------------------------------------------
+ subroutine openunits
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+ character*16 form,nodename
+ integer nodelen
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.CONTROL'
+ integer lenpre,lenpot,ilen,lentmp
+ external ilen
+ character*3 out1file_text,ucase
+ character*3 ll
+ external ucase
+c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
+ call getenv_loc("PREFIX",prefix)
+ pref_orig = prefix
+ call getenv_loc("POT",pot)
+ call getenv_loc("DIRTMP",tmpdir)
+ call getenv_loc("CURDIR",curdir)
+ call getenv_loc("OUT1FILE",out1file_text)
+c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
+ out1file_text=ucase(out1file_text)
+ if (out1file_text(1:1).eq."Y") then
+ out1file=.true.
+ else
+ out1file=fg_rank.gt.0
+ endif
+ lenpre=ilen(prefix)
+ lenpot=ilen(pot)
+ lentmp=ilen(tmpdir)
+ if (lentmp.gt.0) then
+ write (*,'(80(1h!))')
+ write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
+ write (*,'(80(1h!))')
+ write (*,*)"All output files will be on node /tmp directory."
+#ifdef MPI
+ call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
+ if (me.eq.king) then
+ write (*,*) "The master node is ",nodename
+ else if (fg_rank.eq.0) then
+ write (*,*) "I am the CG slave node ",nodename
+ else
+ write (*,*) "I am the FG slave node ",nodename
+ endif
+#endif
+ PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
+ lenpre = lentmp+lenpre+1
+ endif
+ entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
+C Get the names and open the input files
+#if defined(WINIFL) || defined(WINPGI)
+ open(1,file=pref_orig(:ilen(pref_orig))//
+ & '.inp',status='old',readonly,shared)
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+C Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old',readonly,shared)
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old',readonly,shared)
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
+#endif
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old',readonly,shared)
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
+#endif
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old',readonly,shared)
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old',readonly,shared)
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old',readonly,shared)
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old',readonly,shared)
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
+ & action='read')
+c print *,"Processor",myrank," opened file 1"
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+c print *,"Processor",myrank," opened file 9"
+C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+C Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old',action='read')
+c print *,"Processor",myrank," opened file IBOND"
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old',action='read')
+c print *,"Processor",myrank," opened file ITHEP"
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old',action='read')
+#endif
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old',action='read')
+c print *,"Processor",myrank," opened file IROTAM"
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',action='read')
+#endif
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old',action='read')
+c print *,"Processor",myrank," opened file ITORP"
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old',action='read')
+c print *,"Processor",myrank," opened file ITORDP"
+ call getenv_loc('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status='old',action='read')
+c print *,"Processor",myrank," opened file ISCCOR"
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old',action='read')
+c print *,"Processor",myrank," opened file IFOURIER"
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old',action='read')
+c print *,"Processor",myrank," opened file IELEP"
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old',action='read')
+c print *,"Processor",myrank," opened file ISIDEP"
+c print *,"Processor",myrank," opened parameter files"
+#elif (defined G77)
+ open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+C Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old')
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old')
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old')
+#endif
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old')
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old')
+#endif
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old')
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old')
+ call getenv_loc('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status='old')
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old')
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old')
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old')
+#else
+ open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
+ & readonly)
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+C Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old',readonly)
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old',readonly)
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ print *,"thetname_pdb ",thetname_pdb
+ open (ithep_pdb,file=thetname_pdb,status='old',readonly)
+ print *,ithep_pdb," opened"
+#endif
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old',readonly)
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',readonly)
+#endif
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old',readonly)
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old',readonly)
+ call getenv_loc('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status='old',readonly)
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old',readonly)
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old',readonly)
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old',readonly)
+#endif
+#ifndef OLDSCP
+C
+C 8/9/01 In the newest version SCp interaction constants are read from a file
+C Use -DOLDSCP to use hard-coded constants instead.
+C
+ call getenv_loc('SCPPAR',scpname)
+#if defined(WINIFL) || defined(WINPGI)
+ open (iscpp,file=scpname,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open (iscpp,file=scpname,status='old',action='read')
+#elif (defined G77)
+ open (iscpp,file=scpname,status='old')
+#else
+ open (iscpp,file=scpname,status='old',readonly)
+#endif
+#endif
+ call getenv_loc('PATTERN',patname)
+#if defined(WINIFL) || defined(WINPGI)
+ open (icbase,file=patname,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open (icbase,file=patname,status='old',action='read')
+#elif (defined G77)
+ open (icbase,file=patname,status='old')
+#else
+ open (icbase,file=patname,status='old',readonly)
+#endif
+#ifdef MPI
+C Open output file only for CG processes
+c print *,"Processor",myrank," fg_rank",fg_rank
+ if (fg_rank.eq.0) then
+
+ if (nodes.eq.1) then
+ npos=3
+ else
+ npos = dlog10(dfloat(nodes-1))+1
+ endif
+ if (npos.lt.3) npos=3
+ write (liczba,'(i1)') npos
+ form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
+ & //')'
+ write (liczba,form) me
+ outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
+ & liczba(:ilen(liczba))
+ intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
+ & //'.int'
+ pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
+ & //'.pdb'
+ mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
+ & liczba(:ilen(liczba))//'.mol2'
+ statname=prefix(:lenpre)//'_'//pot(:lenpot)//
+ & liczba(:ilen(liczba))//'.stat'
+ if (lentmp.gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
+ & //liczba(:ilen(liczba))//'.stat')
+ rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
+ & //'.rst'
+ if(usampl) then
+ qname=prefix(:lenpre)//'_'//pot(:lenpot)//
+ & liczba(:ilen(liczba))//'.const'
+ endif
+
+ endif
+#else
+ outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
+ intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
+ pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
+ mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
+ statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
+ if (lentmp.gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
+ & //'.stat')
+ rest2name=prefix(:ilen(prefix))//'.rst'
+ if(usampl) then
+ qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
+ endif
+#endif
+#if defined(AIX) || defined(PGI)
+ if (me.eq.king .or. .not. out1file)
+ & open(iout,file=outname,status='unknown')
+#ifdef DEBUG
+ if (fg_rank.gt.0) then
+ write (liczba,'(i3.3)') myrank/nfgtasks
+ write (ll,'(bz,i3.3)') fg_rank
+ open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
+ & status='unknown')
+ endif
+#endif
+ if(me.eq.king) then
+ open(igeom,file=intname,status='unknown',position='append')
+ open(ipdb,file=pdbname,status='unknown')
+ open(imol2,file=mol2name,status='unknown')
+ open(istat,file=statname,status='unknown',position='append')
+ else
+c1out open(iout,file=outname,status='unknown')
+ endif
+#else
+ if (me.eq.king .or. .not.out1file)
+ & open(iout,file=outname,status='unknown')
+#ifdef DEBUG
+ if (fg_rank.gt.0) then
+ write (liczba,'(i3.3)') myrank/nfgtasks
+ write (ll,'(bz,i3.3)') fg_rank
+ open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
+ & status='unknown')
+ endif
+#endif
+ if(me.eq.king) then
+ open(igeom,file=intname,status='unknown',access='append')
+ open(ipdb,file=pdbname,status='unknown')
+ open(imol2,file=mol2name,status='unknown')
+ open(istat,file=statname,status='unknown',access='append')
+ else
+c1out open(iout,file=outname,status='unknown')
+ endif
+#endif
+csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
+csa csa_seed=prefix(:lenpre)//'.CSA.seed'
+csa csa_history=prefix(:lenpre)//'.CSA.history'
+csa csa_bank=prefix(:lenpre)//'.CSA.bank'
+csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
+csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
+csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
+csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
+csa csa_int=prefix(:lenpre)//'.int'
+csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
+csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
+csa csa_in=prefix(:lenpre)//'.CSA.in'
+c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
+C Write file names
+ if (me.eq.king)then
+ write (iout,'(80(1h-))')
+ write (iout,'(30x,a)') "FILE ASSIGNMENT"
+ write (iout,'(80(1h-))')
+ write (iout,*) "Input file : ",
+ & pref_orig(:ilen(pref_orig))//'.inp'
+ write (iout,*) "Output file : ",
+ & outname(:ilen(outname))
+ write (iout,*)
+ write (iout,*) "Sidechain potential file : ",
+ & sidename(:ilen(sidename))
+#ifndef OLDSCP
+ write (iout,*) "SCp potential file : ",
+ & scpname(:ilen(scpname))
+#endif
+ write (iout,*) "Electrostatic potential file : ",
+ & elename(:ilen(elename))
+ write (iout,*) "Cumulant coefficient file : ",
+ & fouriername(:ilen(fouriername))
+ write (iout,*) "Torsional parameter file : ",
+ & torname(:ilen(torname))
+ write (iout,*) "Double torsional parameter file : ",
+ & tordname(:ilen(tordname))
+ write (iout,*) "SCCOR parameter file : ",
+ & sccorname(:ilen(sccorname))
+ write (iout,*) "Bond & inertia constant file : ",
+ & bondname(:ilen(bondname))
+ write (iout,*) "Bending parameter file : ",
+ & thetname(:ilen(thetname))
+ write (iout,*) "Rotamer parameter file : ",
+ & rotname(:ilen(rotname))
+ write (iout,*) "Threading database : ",
+ & patname(:ilen(patname))
+ if (lentmp.ne.0)
+ &write (iout,*)" DIRTMP : ",
+ & tmpdir(:lentmp)
+ write (iout,'(80(1h-))')
+ endif
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine card_concat(card)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ character*(*) card
+ character*80 karta,ucase
+ external ilen
+ read (inp,'(a)') karta
+ karta=ucase(karta)
+ card=' '
+ do while (karta(80:80).eq.'&')
+ card=card(:ilen(card)+1)//karta(:79)
+ read (inp,'(a)') karta
+ karta=ucase(karta)
+ enddo
+ card=card(:ilen(card)+1)//karta
+ return
+ end
+c----------------------------------------------------------------------------------
+ subroutine readrst
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ open(irest2,file=rest2name,status='unknown')
+ read(irest2,*) totT,EK,potE,totE,t_bath
+ do i=1,2*nres
+ read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
+ enddo
+ do i=1,2*nres
+ read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
+ enddo
+ if(usampl) then
+ read (irest2,*) iset
+ endif
+ close(irest2)
+ return
+ end
+c---------------------------------------------------------------------------------
+ subroutine read_fragments
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.CONTROL'
+ read(inp,*) nset,nfrag,npair,nfrag_back
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
+ & " nfrag_back",nfrag_back
+ do iset=1,nset
+ read(inp,*) mset(iset)
+ do i=1,nfrag
+ read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
+ & qinfrag(i,iset)
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
+ & ifrag(2,i,iset), qinfrag(i,iset)
+ enddo
+ do i=1,npair
+ read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
+ & qinpair(i,iset)
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
+ & ipair(2,i,iset), qinpair(i,iset)
+ enddo
+ do i=1,nfrag_back
+ read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
+ & wfrag_back(3,i,iset),
+ & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
+ & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ enddo
+ enddo
+ return
+ end
+c-------------------------------------------------------------------------------
+ subroutine read_dist_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ integer ifrag_(2,100),ipair_(2,100)
+ double precision wfrag_(100),wpair_(100)
+ character*500 controlcard
+c write (iout,*) "Calling read_dist_constr"
+c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+c call flush(iout)
+ call card_concat(controlcard)
+ call readi(controlcard,"NFRAG",nfrag_,0)
+ call readi(controlcard,"NPAIR",npair_,0)
+ call readi(controlcard,"NDIST",ndist_,0)
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+ call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+ call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+ call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+ call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+c write (iout,*) "IFRAG"
+c do i=1,nfrag_
+c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+c enddo
+c write (iout,*) "IPAIR"
+c do i=1,npair_
+c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+c enddo
+ call flush(iout)
+ do i=1,nfrag_
+ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+ if (ifrag_(2,i).gt.nstart_sup+nsup-1)
+ & ifrag_(2,i)=nstart_sup+nsup-1
+c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ call flush(iout)
+ if (wfrag_(i).gt.0.0d0) then
+ do j=ifrag_(1,i),ifrag_(2,i)-1
+ do k=j+1,ifrag_(2,i)
+ write (iout,*) "j",j," k",k
+ ddjk=dist(j,k)
+ if (constr_dist.eq.1) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king)
+ & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,npair_
+ if (wpair_(i).gt.0.0d0) then
+ ii = ipair_(1,i)
+ jj = ipair_(2,i)
+ if (ii.gt.jj) then
+ itemp=ii
+ ii=jj
+ jj=itemp
+ endif
+ do j=ifrag_(1,ii),ifrag_(2,ii)
+ do k=ifrag_(1,jj),ifrag_(2,jj)
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ forcon(nhpb)=wpair_(i)
+ dhpb(nhpb)=dist(j,k)
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king)
+ & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,ndist_
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king)
+ & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+ endif
+ enddo
+ call flush(iout)
+ return
+ end
+c-------------------------------------------------------------------------------
+#ifdef WINIFL
+ subroutine flush(iu)
+ return
+ end
+#endif
+#ifdef AIX
+ subroutine flush(iu)
+ call flush_(iu)
+ return
+ end
+#endif
+c------------------------------------------------------------------------------
+ subroutine copy_to_tmp(source)
+ include "DIMENSIONS"
+ include "COMMON.IOUNITS"
+ character*(*) source
+ character* 256 tmpfile
+ integer ilen
+ external ilen
+ logical ex
+ tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
+ inquire(file=tmpfile,exist=ex)
+ if (ex) then
+ write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
+ & " to temporary directory..."
+ write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
+ call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine move_from_tmp(source)
+ include "DIMENSIONS"
+ include "COMMON.IOUNITS"
+ character*(*) source
+ integer ilen
+ external ilen
+ write (*,*) "Moving ",source(:ilen(source)),
+ & " from temporary directory to working directory"
+ write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
+ call system("/bin/mv "//source(:ilen(source))//" "//curdir)
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine random_init(seed)
+C
+C Initialize random number generator
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef AMD64
+ integer*8 iseedi8
+#endif
+#ifdef MPI
+ include 'mpif.h'
+ logical OKRandom, prng_restart
+ real*8 r1
+ integer iseed_array(4)
+#endif
+ include 'COMMON.IOUNITS'
+ include 'COMMON.TIME1'
+ include 'COMMON.THREAD'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CONTROL'
+ include 'COMMON.MCM'
+ include 'COMMON.MAP'
+ include 'COMMON.HEADER'
+csa include 'COMMON.CSA'
+ include 'COMMON.CHAIN'
+ include 'COMMON.MUCA'
+ include 'COMMON.MD'
+ include 'COMMON.FFIELD'
+ include 'COMMON.SETUP'
+ iseed=-dint(dabs(seed))
+ if (iseed.eq.0) then
+ write (iout,'(/80(1h*)/20x,a/80(1h*))')
+ & 'Random seed undefined. The program will stop.'
+ write (*,'(/80(1h*)/20x,a/80(1h*))')
+ & 'Random seed undefined. The program will stop.'
+#ifdef MPI
+ call mpi_finalize(mpi_comm_world,ierr)
+#endif
+ stop 'Bad random seed.'
+ endif
+#ifdef MPI
+ if (fg_rank.eq.0) then
+ seed=seed*(me+1)+1
+#ifdef AMD64
+ iseedi8=dint(seed)
+ if(me.eq.king .or. .not. out1file)
+ & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
+ write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
+ OKRandom = prng_restart(me,iseedi8)
+#else
+ do i=1,4
+ tmp=65536.0d0**(4-i)
+ iseed_array(i) = dint(seed/tmp)
+ seed=seed-iseed_array(i)*tmp
+ enddo
+ if(me.eq.king .or. .not. out1file)
+ & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
+ & (iseed_array(i),i=1,4)
+ write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
+ & (iseed_array(i),i=1,4)
+ OKRandom = prng_restart(me,iseed_array)
+#endif
+ if (OKRandom) then
+ r1=ran_number(0.0D0,1.0D0)
+ if(me.eq.king .or. .not. out1file)
+ & write (iout,*) 'ran_num',r1
+ if (r1.lt.0.0d0) OKRandom=.false.
+ endif
+ if (.not.OKRandom) then
+ write (iout,*) 'PRNG IS NOT WORKING!!!'
+ print *,'PRNG IS NOT WORKING!!!'
+ if (me.eq.0) then
+ call flush(iout)
+ call mpi_abort(mpi_comm_world,error_msg,ierr)
+ stop
+ else
+ write (iout,*) 'too many processors for parallel prng'
+ write (*,*) 'too many processors for parallel prng'
+ call flush(iout)
+ stop
+ endif
+ endif
+ endif
+#else
+ call vrndst(iseed)
+ write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
+#endif
+ return
+ end
projects / unres.git / commitdiff