+++ /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)
- if (modecalc.eq.8) then
- inquire (file="fort.40",exist=file_exist)
- if (.not.file_exist) call csaread
- 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'
- 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
- 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
- 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(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----------------------------------------------------------------------------
- subroutine csaread
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.IOUNITS'
- include 'COMMON.GEO'
- include 'COMMON.CSA'
- include 'COMMON.BANK'
- include 'COMMON.CONTROL'
- character*80 ucase
- character*620 mcmcard
- call card_concat(mcmcard)
-
- call readi(mcmcard,'NCONF',nconf,50)
- call readi(mcmcard,'NADD',nadd,0)
- call readi(mcmcard,'JSTART',jstart,1)
- call readi(mcmcard,'JEND',jend,1)
- call readi(mcmcard,'NSTMAX',nstmax,500000)
- call readi(mcmcard,'N0',n0,1)
- call readi(mcmcard,'N1',n1,6)
- call readi(mcmcard,'N2',n2,4)
- call readi(mcmcard,'N3',n3,0)
- call readi(mcmcard,'N4',n4,0)
- call readi(mcmcard,'N5',n5,0)
- call readi(mcmcard,'N6',n6,10)
- call readi(mcmcard,'N7',n7,0)
- call readi(mcmcard,'N8',n8,0)
- call readi(mcmcard,'N9',n9,0)
- call readi(mcmcard,'N14',n14,0)
- call readi(mcmcard,'N15',n15,0)
- call readi(mcmcard,'N16',n16,0)
- call readi(mcmcard,'N17',n17,0)
- call readi(mcmcard,'N18',n18,0)
-
- vdisulf=(index(mcmcard,'DYNSS').gt.0)
-
- call readi(mcmcard,'NDIFF',ndiff,2)
- call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
- call readi(mcmcard,'IS1',is1,1)
- call readi(mcmcard,'IS2',is2,8)
- call readi(mcmcard,'NRAN0',nran0,4)
- call readi(mcmcard,'NRAN1',nran1,2)
- call readi(mcmcard,'IRR',irr,1)
- call readi(mcmcard,'NSEED',nseed,20)
- call readi(mcmcard,'NTOTAL',ntotal,10000)
- call reada(mcmcard,'CUT1',cut1,2.0d0)
- call reada(mcmcard,'CUT2',cut2,5.0d0)
- call reada(mcmcard,'ESTOP',estop,-3000.0d0)
- call readi(mcmcard,'ICMAX',icmax,3)
- call readi(mcmcard,'IRESTART',irestart,0)
-c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
- ntbankm=0
-c!bankt
- call reada(mcmcard,'DELE',dele,20.0d0)
- call reada(mcmcard,'DIFCUT',difcut,720.0d0)
- call readi(mcmcard,'IREF',iref,0)
- call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
- call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
- call readi(mcmcard,'NCONF_IN',nconf_in,0)
- call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
- write (iout,*) "NCONF_IN",nconf_in
- return
- 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)
- call getenv_loc('ROTPAR',rotname)
- open (irotam,file=rotname,status='old',readonly,shared)
- 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"
- call getenv_loc('ROTPAR',rotname)
- open (irotam,file=rotname,status='old',action='read')
-c print *,"Processor",myrank," opened file IROTAM"
- 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')
- call getenv_loc('ROTPAR',rotname)
- open (irotam,file=rotname,status='old')
- 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)
- call getenv_loc('ROTPAR',rotname)
- open (irotam,file=rotname,status='old',readonly)
- 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_rbank=prefix(:lenpre)//'.CSA.rbank'
- csa_seed=prefix(:lenpre)//'.CSA.seed'
- csa_history=prefix(:lenpre)//'.CSA.history'
- csa_bank=prefix(:lenpre)//'.CSA.bank'
- csa_bank1=prefix(:lenpre)//'.CSA.bank1'
- csa_alpha=prefix(:lenpre)//'.CSA.alpha'
- csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
-c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
- csa_int=prefix(:lenpre)//'.int'
- csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
- csa_native_int=prefix(:lenpre)//'.CSA.native.int'
- 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'
- 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