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' include 'COMMON.SPLITELE' 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.INTERACT' include 'COMMON.SETUP' include 'COMMON.SPLITELE' include 'COMMON.SHIELD' 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) C this variable with_theta_constr is the variable which allow to read and execute the C constrains on theta angles WITH_THETA_CONSTR is the keyword with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 write (iout,*) "with_theta_constr ",with_theta_constr call readi(controlcard,'SYM',symetr,1) call reada(controlcard,'TIMLIM',timlim,2800.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) AFMlog=(index(controlcard,'AFM')) selfguide=(index(controlcard,'SELFGUIDE')) print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'IPRINT',iprint,0) C SHIELD keyword sets if the shielding effect of side-chains is used C 0 denots no shielding is used all peptide are equally despite the C solvent accesible area C 1 the newly introduced function C 2 reseved for further possible developement call readi(controlcard,'SHIELD',shield_mode,0) C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "shield_mode",shield_mode C endif call readi(controlcard,'TORMODE',tor_mode,0) C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "torsional and valence angle mode",tor_mode 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 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 C DISTCHAINMAX become obsolete for periodic boundry condition call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0) C Reading the dimensions of box in x,y,z coordinates call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) call reada(controlcard,"LIPTHICK",lipthick,0.0d0) call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) if (lipthick.gt.0.0d0) then bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick C endif if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick if ((lipbufthick*2.0d0).gt.lipthick) &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" endif write(iout,*) "bordliptop=",bordliptop write(iout,*) "bordlipbot=",bordlipbot write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot write (iout,*) "SHIELD MODE",shield_mode if (shield_mode.gt.0) then pi=3.141592d0 C VSolvSphere the volume of solving sphere C print *,pi,"pi" C rpp(1,1) is the energy r0 for peptide group contact and will be used for it C there will be no distinction between proline peptide group and normal peptide C group in case of shielding parameters VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 write (iout,*) VSolvSphere,VSolvSphere_div C long axis of side chain do i=1,ntyp long_r_sidechain(i)=vbldsc0(1,i) short_r_sidechain(i)=sigma0(i) enddo buff_shield=1.0d0 endif if (me.eq.king .or. .not.out1file ) & write (iout,*) "DISTCHAINMAX",distchainmax 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 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 write (iout,*) 'Total number of replicas ',iremd_m_total 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 c call reada(controlcard,"R_CUT",r_cut,2.0d0) c call reada(controlcard,"LAMBDA",rlamb,0.3d0) rest = index(controlcard,"REST").gt.0 tbf = index(controlcard,"TBF").gt.0 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(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.SETUP' include 'COMMON.SHIELD' character*4 sequence(maxres) integer rescode double precision x(maxvar) character*256 pdbfile character*400 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. 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) call reada(weightcard,'WSHIELD',wshield,1.0d0) call reada(weightcard,'WLT',wliptran,0.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 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,ntyp 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) call reada(weightcard,"ATRISS",atriss,0.301D0) call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) call reada(weightcard,"DTRISS",dtriss,1.001D0) write (iout,*) "ATRISS=", atriss write (iout,*) "BTRISS=", btriss write (iout,*) "CTRISS=", ctriss write (iout,*) "DTRISS=", dtriss 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 if (wstrain.ne.0.0) then ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain else ss_depth=0.0 endif endif write (iout,*) "wshield,", wshield 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 write (iout,*) 'Begin reading pdb data' c call flush(iout) call readpdb c write (iout,*) 'Finished reading pdb data' c call flush(iout) 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 .and. itype(i).ne.ntyp1) 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(iabs(itype(i))) vbld_inv(i+nres)=dsc_inv(iabs(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.ntyp1 .or. itype(i+1).eq.ntyp1) then #else if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR else if (iabs(itype(i+1)).ne.20) then #else else if (iabs(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 C read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(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,3f8.3)') idih_constr(i),phi0(i),drange(i), & ftors(i) enddo endif do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo C if(me.eq.king.or..not.out1file) C & 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 C first setting the theta boundaries to 0 to pi C this mean that there is no energy penalty for any angle occuring this can be applied C for generate random conformation but is not implemented in this way C do i=1,nres C thetabound(1,i)=0 C thetabound(2,i)=pi C enddo C begin reading theta constrains this is quartic constrains allowing to C have smooth second derivative if (with_theta_constr) then C with_theta_constr is keyword allowing for occurance of theta constrains read (inp,*) ntheta_constr C ntheta_constr is the number of theta constrains if (ntheta_constr.gt.0) then C read (inp,*) ftors read (inp,*) (itheta_constr(i),theta_constr0(i), & theta_drange(i),for_thet_constr(i), & i=1,ntheta_constr) C the above code reads from 1 to ntheta_constr C itheta_constr(i) residue i for which is theta_constr C theta_constr0 the global minimum value C theta_drange is range for which there is no energy penalty C for_thet_constr is the force constant for quartic energy penalty C E=k*x**4 if(me.eq.king.or..not.out1file)then write (iout,*) & 'There are',ntheta_constr,' constraints on phi angles.' do i=1,ntheta_constr write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), & theta_drange(i), & for_thet_constr(i) enddo endif do i=1,ntheta_constr theta_constr0(i)=deg2rad*theta_constr0(i) theta_drange(i)=deg2rad*theta_drange(i) enddo C if(me.eq.king.or..not.out1file) C & write (iout,*) 'FTORS',ftors C do i=1,ntheta_constr C ii = itheta_constr(i) C thetabound(1,ii) = phi0(i)-drange(i) C thetabound(2,ii) = phi0(i)+drange(i) C enddo endif ! ntheta_constr.gt.0 endif! with_theta_constr C C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 C write (iout,*) "with_dihed_constr ",with_dihed_constr 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.ntyp1) nnt=2 if (itype(nres).eq.ntyp1) 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_extconf call setup_var czscore call geom_to_var(nvar,coord_exp_zs(1,1)) nstart_sup=nnt nstart_seq=nnt nsup=nct-nnt+1 kkk=1 do i=1,2*nres do j=1,3 cref(j,i,kkk)=c(j,i) enddo enddo call contact(.true.,ncont_ref,icont_ref,co) endif endif print *, "A TU" write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp 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 C 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) write (iout,*) "Exit READ_CART" write (iout,'(8f10.5)') & ((c(l,k),l=1,3),k=1,nres), & ((c(l,k+nres),l=1,3),k=nnt,nct) call int_from_cart1(.true.) write (iout,*) "Finish INT_TO_CART" 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 .and. itype(i).ne.ntyp1) 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) call chainbuild_extconf 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 if (itype(i).le.0) omeg(i)=-omeg(i) enddo call chainbuild_extconf 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 write (*,'(a)') & ' error in generating random conformation.' stop 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" 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) 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 (ns.gt.0.and.dyn_ss) then 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 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 print *,"A TU?" 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(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ', & restyp(itype(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) if(fg_rank.eq.0) & 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 .and. itype(i).ne.ntyp1) 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 .and. itype(i).ne.ntyp1) 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) 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) 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('TORKCC',torkccname) open (itorkcc,file=torkccname,status='old',readonly,shared) call getenv_loc('THETKCC',thetkccname) open (ithetkcc,file=thetkccname,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) call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,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') call getenv_loc('TORKCC',torkccname) open (itorkcc,file=torkccname,status='old',action='read') call getenv_loc('THETKCC',thetkccname) open (ithetkcc,file=thetkccname,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') call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,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('TORKCC',torkccname) open (itorkcc,file=torkccname,status='old') call getenv_loc('THETKCC',thetkccname) open (ithetkcc,file=thetkccname,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') call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,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('TORKCC',torkccname) open (itorkcc,file=torkccname,status='old',readonly) call getenv_loc('THETKCC',thetkccname) open (ithetkcc,file=thetkccname,status='old',readonly) call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old',readonly) #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) print *,"thetname_pdb ",thetname_pdb open (ithep_pdb,file=thetname_pdb,status='old',action='read') print *,ithep_pdb," opened" #endif 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) call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif #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,*) "Thetpdb parameter file : ", & thetname_pdb(:ilen(thetname_pdb)) 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 totTafm=totT 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_afminp 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' character*320 afmcard print *, "wchodze" call card_concat(afmcard) call readi(afmcard,"BEG",afmbeg,0) call readi(afmcard,"END",afmend,0) call reada(afmcard,"FORCE",forceAFMconst,0.0d0) call reada(afmcard,"VEL",velAFMconst,0.0d0) print *,'FORCE=' ,forceAFMconst CCCC NOW PROPERTIES FOR AFM distafminit=0.0d0 do i=1,3 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit enddo distafminit=dsqrt(distafminit) print *,'initdist',distafminit 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 print *, "WCHODZE" 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) c 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 print *,ndist_ do i=1,ndist_ if (constr_dist.eq.11) then read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) else C print *,"in else" read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1) endif 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 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) C if (forcon(nhpb+1).gt.0.0d0) then C nhpb=nhpb+1 C 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 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 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 if(me.eq.king) & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed OKRandom = prng_restart(me,iseed) #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) & 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 c r1 = prng_next(me) r1=ran_number(0.0D0,1.0D0) if(me.eq.king) & 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