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,i2,3f10.5)') i-nss,ihpb(i),jhpb(i), & ibecarb(i),dhpb(i),dhpb1(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,2) 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) dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. enddo do i=1,maxres-1 do j=i+1,maxres dyn_ssbond_ij(i,j)=1.0d300 enddo enddo call reada(weightcard,"HT",Ht,0.0D0) if (dyn_ss) then ss_depth=ebr/wsc-0.25*eps(1,1) Ht=Ht/wsc-0.25*eps(1,1) akcm=akcm*wstrain/wsc akth=akth*wstrain/wsc akct=akct*wstrain/wsc v1ss=v1ss*wstrain/wsc v2ss=v2ss*wstrain/wsc v3ss=v3ss*wstrain/wsc else ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain endif 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," SS_DEPTH",ss_depth write (iout,*)" HT",Ht 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 C Following 2 lines for diagnostics; comment out if not needed write (iout,*) "Before sideadd" call intout 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 C 10/03/12 Adam: Recalculate coordinates with new side chain positions call chainbuild endif C Following 2 lines for diagnostics; comment out if not needed write (iout,*) "After sideadd" call intout 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 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 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup if (constr_dist.gt.0) then call read_dist_constr call hpb_partition endif c write (iout,*) "After read_dist_constr nhpb",nhpb c call flush(iout) 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) if (dyn_ss) then write(iout,*)"Running with dynamic disulfide-bond formation" do i=nss+1,nhpb ihpb(i-nss)=ihpb(i) jhpb(i-nss)=jhpb(i) forcon(i-nss)=forcon(i) dhpb(i-nss)=dhpb(i) enddo nhpb=nhpb-nss nss=0 call hpb_partition do i=1,ns dyn_ss_mask(iss(i))=.true. enddo else 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 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) print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1) print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1) print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1) 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', &action='read') 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',action='read') call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) print *,"thetname_pdb ",thetname_pdb open (ithep_pdb,file=thetname_pdb,status='old',action='read') print *,ithep_pdb," opened" #endif call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',action='read') #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') call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,status='old',action='read') call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old',action='read') call getenv_loc('FOURIER',fouriername) open (ifourier,file=fouriername,status='old',action='read') call getenv_loc('ELEPAR',elename) open (ielep,file=elename,status='old',action='read') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',action='read') #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',action='read') #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',action='read') #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') c#define DEBUG #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 c#undef DEBUG 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') c#define DEBUG #ifdef DEBUG if (fg_rank.gt.0) then print "Processor",fg_rank," opening output file" 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 c#undef DEBUG 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 if (.not.refstr .and. nfrag.gt.0) then write (iout,*) & "ERROR: no reference structure to compute distance restraints" write (iout,*) & "Restraints must be specified explicitly (NDIST=number)" stop endif if (nfrag.lt.2 .and. npair.gt.0) then write (iout,*) "ERROR: Less than 2 fragments specified", & " but distance restraints between pairs requested" stop endif 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),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 if (ibecarb(i).gt.0) then ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres endif if (dhpb(nhpb).eq.0.0d0) & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) endif enddo #ifdef MPI if (.not.out1file .or. me.eq.king) then #endif do i=1,nhpb write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) enddo call flush(iout) #ifdef MPI endif #endif 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