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