From: Adam Liwo Date: Thu, 26 Nov 2015 21:05:20 +0000 (+0100) Subject: Fixed eello5, eello6, eturn6, and shortrange RESPA X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=e4035e50115ab9c0e65b0445aed96096a5cb86d5 Fixed eello5, eello6, eturn6, and shortrange RESPA --- diff --git a/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe b/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe index 931588b..9d96d80 100755 Binary files a/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe and b/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe differ diff --git a/source/unres/src_MD-M/Makefile b/source/unres/src_MD-M/Makefile deleted file mode 100644 index 35c2a1f..0000000 --- a/source/unres/src_MD-M/Makefile +++ /dev/null @@ -1,151 +0,0 @@ -#INSTALL_DIR = /usr/local/mpich-1.2.0 -INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh -# -#FC= /usr/local/opt/intel/compiler60/ia32/bin/ifc -FC= ifort - -OPT = -O3 -ip -w -#OPT = -g -CB -#OPT = -g -CFLAGS = -DSGI -c - -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-M/unres_Tc_procor_newparm_em64-D-symetr.exe -#LIBS = -L$(INSTALL_DIR)/lib_pgi -lmpich xdrf/libxdrf.a -#LIBS = -L$(INSTALL_DIR)/lib_ifort -lmpich xdrf/libxdrf.a -LIBS = -L$(INSTALL_DIR)/lib -lmpich ../../lib/xdrf_em64/libxdrf.a - -ARCH = LINUX -PP = /lib/cpp -P - - -all: GAB - -.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_CSA.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 \ - MP.o compare_s1.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 rattle.o gauss.o energy_split-sep.o \ - q_measure.o gnmr1.o test.o ssMD.o permut.o distfit.o checkvar.o - -no_option: - -GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI -DAMD64 \ - -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC\ - -DSCCORPDB -GAB: BIN = ../../../bin/unres/MD-M/unres_ifort_MPICH_GAB.exe -GAB: ${object} ../../lib/xdrf_em64/libxdrf.a - cc -o compinfo compinfo.c - ./compinfo | true - ${FC} ${FFLAGS} cinfo.f - ${FC} ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} - -E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI -DAMD64 \ - -DSPLITELE -DLANG0 -E0LL2Y: BIN = ../../../bin/unres/MD-M/unres_ifort_MPICH_E0LL2Y.exe -E0LL2Y: ${object} ../../lib/xdrf_em64/libxdrf.a - cc -o compinfo compinfo.c - ./compinfo | true - ${FC} ${FFLAGS} cinfo.f - ${FC} ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} - -../../lib/xdrf_em64/libxdrf.a: - cd ../../lib/xdrf_em64 && make - -4P: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN \ - -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -4P: BIN = ../../../bin/unres/MD/unres-mult_ifort_single_4P.exe -4P: ${object} xdrf/libxdrf.a - cc -o compinfo compinfo.c - ./compinfo | true - ${FC} ${FFLAGS} cinfo.f - ${FC} ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} - -E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN \ - -DSPLITELE -DLANG0 -E0LL2Y: BIN = ../../../bin/unres/MD/unres-mult_ifort_single_E0LL2Y.exe -E0LL2Y: ${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 -f *.o && /bin/rm -f compinfo && cd xdrf && make clean - -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-M/Makefile b/source/unres/src_MD-M/Makefile new file mode 120000 index 0000000..8453cdd --- /dev/null +++ b/source/unres/src_MD-M/Makefile @@ -0,0 +1 @@ +Makefile_MPICH_ifort \ No newline at end of file diff --git a/source/unres/src_MD-M/cinfo.f b/source/unres/src_MD-M/cinfo.f index c7b5141..8757613 100644 --- a/source/unres/src_MD-M/cinfo.f +++ b/source/unres/src_MD-M/cinfo.f @@ -1,15 +1,15 @@ C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C -C 0 40360 9 +C 0 40360 15 subroutine cinfo include 'COMMON.IOUNITS' write(iout,*)'++++ Compile info ++++' - write(iout,*)'Version 0.40360 build 9' - write(iout,*)'compiled Fri Jan 23 21:00:08 2015' - write(iout,*)'compiled by adam@mmka' + write(iout,*)'Version 0.40360 build 15' + write(iout,*)'compiled Thu Nov 26 21:54:36 2015' + write(iout,*)'compiled by adam@piasek4' write(iout,*)'OS name: Linux ' - write(iout,*)'OS release: 3.2.0-72-generic ' + write(iout,*)'OS release: 3.2.0-91-generic ' write(iout,*)'OS version:', - & ' #107-Ubuntu SMP Thu Nov 6 14:24:01 UTC 2014 ' + & ' #129-Ubuntu SMP Wed Sep 9 10:56:06 UTC 2015 ' write(iout,*)'flags:' write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...' write(iout,*)'FC= ifort' diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index c0b4c84..1c6470d 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -885,14 +885,16 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) +c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +c call flush(iout) ind=ind+1 itypj=itype(j) if (itypj.eq.ntyp1) cycle c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) -c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, +c write (iout,*) "j",j,itypi,itypj,dsc_inv(itypj),dscj_inv, c & 1.0d0/vbld(j+nres) -c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +c call flush(iout) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@ -911,6 +913,8 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) +c write (iout,*) "After box" +c call flush(iout) if (zj.lt.0) zj=zj+boxzsize if ((zj.gt.bordlipbot) &.and.(zj.lt.bordliptop)) then @@ -934,6 +938,8 @@ C lipbufthick is thickenes of lipid buffore sslipj=0.0d0 ssgradlipj=0.0 endif +c write (iout,*) "After lipid" +c call flush(iout) aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 @@ -981,6 +987,8 @@ C lipbufthick is thickenes of lipid buffore C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular +c write (iout,*) "After sc_angular" +c call flush(iout) sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij @@ -989,9 +997,9 @@ c rij_shift=1.2*sig0ij C I hate to put IF's in the loops, but here don't have another choice!!!! if (rij_shift.le.0.0D0) then evdw=1.0D20 -cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') -cd & restyp(itypi),i,restyp(itypj),j, -cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) +c write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c & restyp(itypi),i,restyp(itypj),j, +c & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) return endif sigder=-sig*sigsq @@ -1005,6 +1013,7 @@ c--------------------------------------------------------------- eps3der=evdwij*eps2rt c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 +c call flush(iout) evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij*sss if (lprn) then @@ -1016,10 +1025,13 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, & evdwij + call flush(iout) endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij + if (energy_dec) then + write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij + call flush(iout) + endif C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1034,8 +1046,12 @@ C Calculate the radial part of the gradient gg(3)=zj*fac gg_lipi(3)=ssgradlipi*evdwij gg_lipj(3)=ssgradlipj*evdwij +c write (iout,*) "Calling sc_grad_scale" +c call flush(iout) C Calculate angular part of the gradient. call sc_grad_scale(sss) +c write (iout,*) "After sc_grad_scale" +c call flush(iout) endif enddo ! j enddo ! iint @@ -1298,6 +1314,8 @@ C---------------------------------------------------------------------------- include 'COMMON.IOUNITS' double precision dcosom1(3),dcosom2(3) double precision scalfac +c write (iout,*) "sc_grad_scale",i,j,k,l +c call flush(iout) eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 @@ -1335,8 +1353,8 @@ C C Calculate the components of the gradient in DC and X C do l=1,3 - gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k) - gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k) + gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) enddo return end diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 00862f2..747b132 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -130,6 +130,8 @@ cmc c if (dyn_ss) call dyn_set_nss c print *,"Processor",myrank," computed USCSC" +c write (iout,*) "SCSC computed OK" +c call flush_(iout) #ifdef TIMING time01=MPI_Wtime() #endif @@ -171,6 +173,8 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, c & eello_turn4) endif +c write (iout,*) "eelec computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed UELEC" C C Calculate excluded-volume interaction energy between peptide groups @@ -187,10 +191,14 @@ C c write (iout,*) "Soft-sphere SCP potential" call escp_soft_sphere(evdw2,evdw2_14) endif +c write (iout,*) "escp computed OK" +c call flush_(iout) c c Calculate the bond-stretching energy c call ebond(estr) +c write (iout,*) "ebond computed OK" +c call flush_(iout) C C Calculate the disulfide-bridge and other energy and the contributions C from other distance constraints. @@ -206,12 +214,16 @@ C ebe=0 ethetacnstr=0 endif +c write (iout,*) "ebend computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed UB" C C Calculate the SC local energy. C C print *,"TU DOCHODZE?" call esc(escloc) +c write (iout,*) "esc computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed USC" C C Calculate the virtual-bond torsional energy. @@ -223,6 +235,8 @@ cd print *,'nterm=',nterm etors=0 edihcnstr=0 endif +c write (iout,*) "etor computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Utor" C C 6/23/01 Calculate double-torsional energy @@ -232,6 +246,8 @@ C else etors_d=0 endif +c write (iout,*) "etor_d computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Utord" C C 21/5/07 Calculate local sicdechain correlation energy @@ -241,6 +257,8 @@ C else esccor=0.0d0 endif +c write (iout,*) "eback_sc_corr computed OK" +c call flush_(iout) C print *,"PRZED MULIt" c print *,"Processor",myrank," computed Usccorr" C @@ -250,9 +268,13 @@ C n_corr1=0 if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 & .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then +c write (iout,*) "Calling multibody_eello" +c call flush_(iout) call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1) -cd write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1, -cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 +c write(iout,*) +c & 'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1, +c & " ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 +c call flush_(iout) else ecorr=0.0d0 ecorr5=0.0d0 @@ -260,9 +282,14 @@ cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 eturn6=0.0d0 endif if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then +c write (iout,*) "Calling multibody_gb_ecorr" +c call flush_(iout) call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) -cd write (iout,*) "multibody_hb ecorr",ecorr +c write (iout,*) "Exited multibody_hb ecorr",ecorr +c call flush_(iout) endif +c write (iout,*) "multibody computed OK" +c call flush_(iout) c print *,"Processor",myrank," computed Ucorr" C C If performing constraint dynamics, call the constraint energy @@ -281,12 +308,15 @@ C print *,"przed lipidami" if (wliptran.gt.0) then call Eliptransfer(eliptran) endif -C print *,"za lipidami" +c write (iout,*) "lipid energy computed OK" +c call flush_(iout) if (AFMlog.gt.0) then call AFMforce(Eafmforce) else if (selfguide.gt.0) then call AFMvel(Eafmforce) endif +c write (iout,*) "AFMforce computed OK" +c call flush_(iout) #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -335,7 +365,11 @@ c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) +c write (iout,*) "sum energy OK" +c call flush_(iout) if (dyn_ss) call dyn_set_nss +c write (iout,*) "Exiting energy" +c call flush_(iout) c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 @@ -5638,6 +5672,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end +c write (iout,*) "ebend: i=",i +c call flush_(iout) if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. @@ -5739,9 +5775,12 @@ C Derivatives of the "mean" values in gamma1 and gamma2. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo +c write (iout,*) "Exit loop" +c call flush_(iout) ethetacnstr=0.0d0 -C print *,ithetaconstr_start,ithetaconstr_end,"TU" - do i=ithetaconstr_start,ithetaconstr_end +c write (iout,*) ithetaconstr_start,ithetaconstr_end,"TU" +c call flush_(iout) + do i=max0(ithetaconstr_start,1),ithetaconstr_end itheta=itheta_constr(i) thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) @@ -5766,6 +5805,8 @@ C print *,ithetaconstr_start,ithetaconstr_end,"TU" & gloc(itheta+nphi-2,icg) endif enddo +c write (iout,*) "Exit ebend" +c call flush_(iout) C Ufff.... We've done all this!!! return @@ -6103,7 +6144,7 @@ c lprn1=.false. C now constrains ethetacnstr=0.0d0 C print *,ithetaconstr_start,ithetaconstr_end,"TU" - do i=ithetaconstr_start,ithetaconstr_end + do i=max0(ithetaconstr_start,1),ithetaconstr_end itheta=itheta_constr(i) thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) @@ -6986,7 +7027,7 @@ c---------------------------------------------------------------------------- logical lprn C Set lprn=.true. for debugging lprn=.false. -c lprn=.true. +c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end C ANY TWO ARE DUMMY ATOMS in row CYCLE @@ -7379,7 +7420,7 @@ C Set lprn=.true. for debugging if (lprn) then write (iout,'(a)') 'Contact function values before RECEIVE:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,f5.2))') + write (iout,'(2i3,50(1x,i3,f5.2))') & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), & j=1,num_cont_hb(i)) enddo @@ -7618,6 +7659,7 @@ C Calculate the local-electrostatic correlation terms jp1=iabs(j1) c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, c & ' jj=',jj,' kk=',kk +c call flush(iout) if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 & .or. j.lt.0 .and. j1.gt.0) .and. & (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then @@ -7635,8 +7677,9 @@ c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0) enddo ! kk do kk=1,num_conti j1=jcont_hb(kk,i) -c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, -c & ' jj=',jj,' kk=',kk +c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, +c & ' jj=',jj,' kk=',kk +c call flush(iout) if (j1.eq.j+1) then C Contacts I-J and (I+1)-J occur simultaneously. C The system loses extra energy. @@ -7730,6 +7773,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.CONTACTS' include 'COMMON.CHAIN' include 'COMMON.CONTROL' + include 'COMMON.TORSION' double precision gx(3),gx1(3) integer num_cont_hb_old(maxres) logical lprn,ldone @@ -7738,6 +7782,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding C Set lprn=.true. for debugging lprn=.false. eturn6=0.0d0 +c write (iout,*) "MULTIBODY_EELLO" +c call flush(iout) #ifdef MPI do i=1,nres num_cont_hb_old(i)=num_cont_hb(i) @@ -7748,12 +7794,12 @@ C Set lprn=.true. for debugging if (lprn) then write (iout,'(a)') 'Contact function values before RECEIVE:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,f5.2))') + write (iout,'(2i3,50(1x,i3,f5.2))') & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), & j=1,num_cont_hb(i)) enddo + call flush(iout) endif - call flush(iout) do i=1,ntask_cont_from ncont_recv(i)=0 enddo @@ -7955,10 +8001,15 @@ c call flush(iout) if (lprn) then write (iout,'(a)') 'Contact function values:' do i=nnt,nct-2 - write (iout,'(2i3,50(1x,i2,5f6.3))') + write (iout,'(2i3,50(1x,i3,5f6.3))') & i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i), & ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i)) enddo + write (iout,*) "itortyp" + do i=1,nres + write (iout,*) i,itype(i),itortyp(itype(i)) + enddo + call flush(iout) endif ecorr=0.0D0 ecorr5=0.0d0 @@ -7999,8 +8050,11 @@ c write (iout,*) "corr loop i",i do kk=1,num_conti1 j1=jcont_hb(kk,i1) jp1=iabs(j1) -c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, -c & ' jj=',jj,' kk=',kk + if (lprn) then + write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, + & ' jj=',jj,' kk=',kk + call flush(iout) + endif c if (j1.eq.j+1 .or. j1.eq.j-1) then if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 & .or. j.lt.0 .and. j1.gt.0) .and. @@ -8040,6 +8094,12 @@ c do iii=1,nres c write (iout,'(i5,3f10.5)') c & iii,(gradcorr5(jjj,iii),jjj=1,3) c enddo +c write (iout,*) "ecorr4" +c call flush(iout) +c write (iout,*) "eello5:",i,jp,i+1,jp1,jj,kk, +c & itype(jp),itype(i+1),itype(jp1), +c & itortyp(itype(jp)),itortyp(itype(i+1)),itortyp(itype(jp1)) +c call flush(iout) if (wcorr5.gt.0.0d0) & ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk) c write (iout,*) "gradcorr5 after eello5" @@ -8052,12 +8112,16 @@ c enddo 2 'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk) cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6 cd write(2,*)'ijkl',i,jp,i+1,jp1 +c write (iout,*) "ecorr5" +c call flush(iout) if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3 & .or. wturn6.eq.0.0d0))then cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1 ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk) if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') 1 'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk) +c write (iout,*) "ecorr6" +c call flush(iout) cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5, cd & 'ecorr6=',ecorr6 cd write (iout,'(4e15.5)') sred_geom, @@ -8071,10 +8135,13 @@ cd write (iout,*) '******eturn6: i,j,i+1,j1',i,jip,i+1,jp1 if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') 1 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk) cd write (2,*) 'multibody_eello:eturn6',eturn6 +c write (iout,*) "ecorr4" +c call flush(iout) endif ENDIF 1111 continue endif + if (energy_dec) call flush(iout) enddo ! kk enddo ! jj enddo ! i @@ -8837,9 +8904,9 @@ cd endif cd write (iout,*) cd & 'EELLO5: Contacts have occurred for peptide groups',i,j, cd & ' and',k,l - itk=itortyp(itype(k)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) +c itk=itortyp(itype(k)) +c itl=itortyp(itype(l)) +c itj=itortyp(itype(j)) eello5_1=0.0d0 eello5_2=0.0d0 eello5_3=0.0d0 @@ -8908,7 +8975,7 @@ C Cartesian gradient c goto 1112 c1111 continue C Contribution from graph II - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8989,7 +9056,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV cd1110 continue - call transpose2(EE(1,1,itl),auxmat(1,1)) + call transpose2(EE(1,1,l),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9062,7 +9129,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV 1110 continue - call transpose2(EE(1,1,itj),auxmat(1,1)) + call transpose2(EE(1,1,j),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9359,7 +9426,7 @@ C o o o o C C i i C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - itk=itortyp(itype(k)) +c itk=itortyp(itype(k)) s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i)) s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k)) s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k)) @@ -9649,19 +9716,19 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. - iti=itortyp(itype(i)) - if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) - else - itj1=ntortyp - endif - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) - if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) - else - itl1=ntortyp - endif +c iti=itortyp(itype(i)) +c if (j.lt.nres-1) then +c itj1=itortyp(itype(j+1)) +c else +c itj1=ntortyp +c endif +c itk=itortyp(itype(k)) +c itk1=itortyp(itype(k+1)) +c if (l.lt.nres-1) then +c itl1=itortyp(itype(l+1)) +c else +c itl1=ntortyp +c endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) #endif @@ -9669,7 +9736,7 @@ C energy moment and not to the cluster cumulant. s2=0.5d0*scalar2(b1(1,k),auxvec(1)) call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1)) s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9767,25 +9834,25 @@ C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. cd write (2,*) 'eello_graph4: wturn6',wturn6 - iti=itortyp(itype(i)) - itj=itortyp(itype(j)) - if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) - else - itj1=ntortyp - endif - itk=itortyp(itype(k)) - if (k.lt.nres-1) then - itk1=itortyp(itype(k+1)) - else - itk1=ntortyp - endif - itl=itortyp(itype(l)) - if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) - else - itl1=ntortyp - endif +c iti=itortyp(itype(i)) +c itj=itortyp(itype(j)) +c if (j.lt.nres-1) then +c itj1=itortyp(itype(j+1)) +c else +c itj1=ntortyp +c endif +c itk=itortyp(itype(k)) +c if (k.lt.nres-1) then +c itk1=itortyp(itype(k+1)) +c else +c itk1=ntortyp +c endif +c itl=itortyp(itype(l)) +c if (l.lt.nres-1) then +c itl1=itortyp(itype(l+1)) +c else +c itl1=ntortyp +c endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, cd & ' itl',itl,' itl1',itl1 @@ -10798,7 +10865,7 @@ C gradient po costhet &*VofOverlap C grad_shield_side is Cbeta sidechain gradient grad_shield_side(j,ishield_list(i),i)= - & (sh_frac_dist_grad(j)*-2.0d0 + & (sh_frac_dist_grad(j)*(-2.0d0) & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet) & +scale_fac_dist*(cosphi_grad_long(j)) & *2.0d0/(1.0-cosphi)) diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F index 24ab8dd..b9fc9e0 100644 --- a/source/unres/src_MD-M/energy_split-sep.F +++ b/source/unres/src_MD-M/energy_split-sep.F @@ -36,7 +36,7 @@ c if (fg_rank.eq.0) call int_from_cart1(.false.) #ifdef MPI c write(iout,*) "ETOTAL_LONG Processor",fg_rank, c & " absolute rank",myrank," nfgtasks",nfgtasks - call flush(iout) +c call flush(iout) if (nfgtasks.gt.1) then time00=MPI_Wtime() C FG slaves call the following matching MPI_Bcast in ERGASTULUM @@ -369,6 +369,8 @@ c call int_from_cart1(.false.) C C Compute the side-chain and electrostatic interaction energy C +c write (iout,*) "Computing SC" +c call flush(iout) goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj_short(evdw) @@ -396,43 +398,66 @@ C c c Calculate the short-range part of Evdwpp c +c write (iout,*) "SC_short computed OK" +c call flush(iout) call evdwpp_short(evdw1) +c write (iout,*) "VDWPP_short computed OK" +c call flush(iout) c c Calculate the short-range part of ESCp c if (ipot.lt.6) then call escp_short(evdw2,evdw2_14) +c write (iout,*) "SCP_short computed OK" +c call flush(iout) endif c c Calculate the bond-stretching energy c call ebond(estr) +c write (iout,*) "ebond_short computed OK" +c call flush(iout) C C Calculate the disulfide-bridge and other energy and the contributions C from other distance constraints. call edis(ehpb) +c write (iout,*) "ehpb_short computed OK" +c call flush(iout) C C Calculate the virtual-bond-angle energy. C - call ebend(ebe) + ethetacnstr=0.0d0 + call ebend(ebe,ethetacnstr) +c write (iout,*) "ebe_short computed OK" +c call flush(iout) C C Calculate the SC local energy. C call vec_and_deriv +c write (iout,*) "vec_and_dervi_short computed OK" +c call flush(iout) call esc(escloc) +c write (iout,*) "SCLOCAL_short computed OK" +c call flush(iout) C C Calculate the virtual-bond torsional energy. C call etor(etors,edihcnstr) +c write (iout,*) "TOR_short computed OK" +c call flush(iout) C C 6/23/01 Calculate double-torsional energy C call etor_d(etors_d) +c write (iout,*) "DTOR_short computed OK" +c call flush(iout) C C 21/5/07 Calculate local sicdechain correlation energy C if (wsccor.gt.0.0d0) then call eback_sc_corr(esccor) +c write (iout,*) "SCCORR_short computed OK" +c call flush(iout) else esccor=0.0d0 endif @@ -464,9 +489,9 @@ C energia(19)=edihcnstr energia(21)=esccor c write (iout,*) "ETOTAL_SHORT before SUM_ENERGY" - call flush(iout) +c call flush(iout) call sum_energy(energia,.true.) c write (iout,*) "Exit ETOTAL_SHORT" - call flush(iout) +c call flush(iout) return end diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index ae4d710..23fe7df 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -58,7 +58,7 @@ C Assign virtual-bond length vblinv2=vblinv*vblinv c c Read the virtual-bond parameters, masses, and moments of inertia -c and Stokes' radii of the peptide group and side chains +c and Stokes radii of the peptide group and side chains c #ifdef CRYST_BOND read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok @@ -99,13 +99,30 @@ c & vbldsc0(j,i),aksc(j,i),abond0(j,i) enddo enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif C reading lipid parameters - read(iliptranpar,*) pepliptran + read(iliptranpar,*,end=120,err=120) pepliptran do i=1,ntyp - read(iliptranpar,*) liptranene(i) + read(iliptranpar,*,end=120,err=120) liptranene(i) enddo close(iliptranpar) + if (lprint) then + write (iout,*) "Lipid transfer parameters" + write (iout,'(a5,f10.5)') "pept",pepliptran + do i=1,ntyp + write (iout,'(a5,f10.5)') restyp(i),liptranene(i) + enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif + endif #ifdef CRYST_THETA C C Read the parameters of the probability distribution/energy expression @@ -182,6 +199,7 @@ C write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i), & sig0(i),(gthet(j,i),j=1,3) enddo + call flush(iout) else write (iout,'(a)') & 'Parameters of the virtual-bond valence angles:' @@ -210,6 +228,11 @@ C write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i), & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0 enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif endif #else @@ -395,9 +418,13 @@ C enddo enddo enddo - call flush(iout) +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif - write (2,*) "Start reading THETA_PDB",ithep_pdb + write (iout,*) "Start reading THETA_PDB",ithep_pdb do i=1,ntyp c write (2,*) 'i=',i read (ithep_pdb,*,err=111,end=111) @@ -442,7 +469,7 @@ c write (2,*) 'i=',i gthet(j,i)=gthet(j,-i) enddo enddo - write (2,*) "End reading THETA_PDB" + write (iout,*) "End reading THETA_PDB" close (ithep_pdb) #endif close(ithep) @@ -530,6 +557,11 @@ C BSC is amplitude of Gaussian endif endif enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif #else C @@ -550,7 +582,7 @@ C C Read the parameters of the probability distribution/energy expression C of the side chains. C - write (2,*) "Start reading ROTAM_PDB" + write (iout,*) "Start reading ROTAM_PDB" do i=1,ntyp read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i) if (i.eq.10) then @@ -589,7 +621,12 @@ C endif enddo close (irotam_pdb) - write (2,*) "End reading ROTAM_PDB" +c write (iout,*) "End reading ROTAM_PDB" +c#ifdef AIX +c call flush_(iout) +c#else +c call flush(iout) +c#endif #endif close(irotam) @@ -617,6 +654,11 @@ C write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old) enddo enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif #else C @@ -659,8 +701,16 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo enddo close (itorp) +c write (iout,*) "End reading torsional parameters" +c#ifdef AIX +c call flush_(iout) +c#else +c call flush(iout) +c#endif if (lprint) then write (iout,'(/a/)') 'Torsional constants:' + do iblock=1,2 + write (iout,*) "IBLOCK",iblock do i=1,ntortyp do j=1,ntortyp write (iout,*) 'ityp',i,' jtyp',j @@ -676,6 +726,12 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo enddo enddo + enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif C @@ -773,6 +829,11 @@ C Martix of D parameters for two dimesional fourier series enddo enddo enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif #endif C Read of Side-chain backbone correlation parameters @@ -780,6 +841,12 @@ C Modified 11 May 2012 by Adasko CCC C read (isccor,*,end=119,err=119) nsccortyp +c write (iout,*) "Reading sccor parameters",nsccortyp +c#ifdef AIX +c call flush_(iout) +c#else +c call flush(iout) +c#endif #ifdef SCCORPDB read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp) do i=-ntyp,-1 @@ -882,15 +949,23 @@ cc maxinter is maximum interaction sites v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &(1+vlor3sccor(k,i,j)**2) enddo - v0sccor(i,j,iblock)=v0ijsccor + v0sccor(l,i,j)=v0ijsccor enddo enddo enddo close (isccor) #endif +c write (iout,*) "sccor parameters read" +c#ifdef AIX +c call flush_(iout) +c#else +c call flush(iout) +c#endif if (lprint) then - write (iout,'(/a/)') 'Torsional constants:' + write (iout,'(/a/)') 'SC-torsional constants:' + do l=1,maxinter + write (iout,*) "Torsional type",l do i=1,nsccortyp do j=1,nsccortyp write (iout,*) 'ityp',i,' jtyp',j @@ -905,6 +980,12 @@ cc maxinter is maximum interaction sites enddo enddo enddo + enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif C @@ -1047,6 +1128,11 @@ c lprint=.true. write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) enddo enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif c lprint=.false. @@ -1077,6 +1163,11 @@ c lprint=.true. c lprint=.false. enddo enddo +#ifdef AIX + if (lprint) call flush_(iout) +#else + if (lprint) call flush(iout) +#endif C C Read side-chain interaction parameters. C @@ -1251,6 +1342,11 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) endif enddo enddo +#ifdef AIX + if (lprint) call flush_(iout) +#else + if (lprint) call flush(iout) +#endif #ifdef OLDSCP C C Define the SC-p interaction constants (hard-coded; old style) @@ -1297,6 +1393,11 @@ c lprint=.true. write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1), & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2) enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif endif c lprint=.false. #endif @@ -1357,25 +1458,82 @@ C buff_shield=1.0d0 C endif return 111 write (iout,*) "Error reading bending energy parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 112 write (iout,*) "Error reading rotamer energy parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 113 write (iout,*) "Error reading torsional energy parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 114 write (iout,*) "Error reading double torsional energy parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 115 write (iout,*) & "Error reading cumulant (multibody energy) parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 116 write (iout,*) "Error reading electrostatic energy parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip" +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 117 write (iout,*) "Error reading side chain interaction parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 118 write (iout,*) "Error reading SCp interaction parameters." +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif goto 999 119 write (iout,*) "Error reading SCCOR parameters" +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif + goto 999 + 120 write (iout,*) "Error reading lipid parameters" +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif 999 continue #ifdef MPI call MPI_Finalize(Ierror)