From: Adam Liwo Date: Sun, 3 Jan 2016 20:22:38 +0000 (+0100) Subject: New valence-torsionals completed X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=00ff2d632b212c4d4a388e8f7f5394763b65e3bb New valence-torsionals completed Introduced itype2loc table to denote local fourier types; now different from itortyp Fixed gradients near gamma=180 --- diff --git a/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe b/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe deleted file mode 100755 index 25bab3f..0000000 Binary files a/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe and /dev/null differ diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index d03828e..e13259d 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -5,7 +5,9 @@ & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, & gnorm_check,gradout,split_ene,with_theta_constr - common /cntrl/ modecalc,iscode,indpdb,indback,indphi,iranconf, + double precision aincr + common /cntrl/ aincr,modecalc,iscode,indpdb,indback,indphi, + & iranconf, & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint, & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, diff --git a/source/unres/src_MD-M/COMMON.IOUNITS b/source/unres/src_MD-M/COMMON.IOUNITS index 8c4fdb7..28a14a1 100644 --- a/source/unres/src_MD-M/COMMON.IOUNITS +++ b/source/unres/src_MD-M/COMMON.IOUNITS @@ -11,12 +11,13 @@ C General I/O units & files integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam, & itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,istat, & ientin,ientout,izs1,isecpred,ibond,irest2,iifrag,icart, - & irest1,isccor,ithep_pdb,irotam_pdb,itorkcc,ithetkcc + & irest1,isccor,ithep_pdb,irotam_pdb,itorkcc,ithetkcc, + & iliptranpar common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep, & irotam,itorp,itordp,ifourier,ielep,isidep,iscpp,icbase, & istat,ientin,ientout,izs1,isecpred,ibond,irest2,iifrag, & icart,irest1,isccor,ithep_pdb,irotam_pdb, - & itorkcc,ithetkcc + & itorkcc,ithetkcc,iliptranpar character*256 outname,intname,pdbname,mol2name,statname,intinname, & entname,prefix,secpred,rest2name,qname,cartname,tmpdir, & mremd_rst_name,curdir,pref_orig diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION index 56c86ac..e7c54a0 100644 --- a/source/unres/src_MD-M/COMMON.TORSION +++ b/source/unres/src_MD-M/COMMON.TORSION @@ -1,6 +1,6 @@ C Torsional constants of the rotation about virtual-bond dihedral angles double precision v1,v2,vlor1,vlor2,vlor3,v0,v1_kcc,v2_kcc, - & v1_chyb,v2_chyb,v1bend_chyb + & v11_chyb,v21_chyb,v12_chyb,v22_chyb,v1bend_chyb integer itortyp,ntortyp,nterm,nlor,nterm_old,nterm_kcc_Tb, & nterm_kcc,itortyp_kcc,nbend_kcc_Tb common/torsion/v0(-maxtor:maxtor,-maxtor:maxtor,2), @@ -10,8 +10,10 @@ C Torsional constants of the rotation about virtual-bond dihedral angles & vlor2(maxlor,maxtor,maxtor),vlor3(maxlor,maxtor,maxtor), & v1_kcc(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), & v2_kcc(maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), - & v1_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), - & v2_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), + & v11_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), + & v21_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), + & v12_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), + & v22_chyb(maxtermkcc,maxtermkcc,-maxtor:maxtor,-maxtor:maxtor), & v1bend_chyb(maxtermkcc,-maxtor:maxtor), & itortyp(-ntyp1:ntyp1),ntortyp, & itortyp_kcc(-ntyp1:ntyp1), @@ -36,17 +38,16 @@ C 6/23/01 - constants for double torsionals C 9/18/99 - added Fourier coeffficients of the expansion of local energy C surface double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde, - &bnew1,bnew2,eenew,gtb1,gtb2,eeold,gtee - integer nloctyp - common/fourier/ b1(2,0:maxres),b2(2,0:maxres),b(13,0:maxtor), - & bnew1(3,2,-maxtor:maxtor),bnew2(3,2,-maxtor:maxtor), - & cc(2,2,-maxtor:maxtor), - & dd(2,2,-maxtor:maxtor),eeold(2,2,-maxtor:maxtor), - & eenew(2,-maxtor:maxtor), + & b,bnew1,bnew2,eenew,gtb1,gtb2,eeold,gtee + integer nloctyp,iloctyp(-ntyp1:ntyp1),itype2loc(-ntyp1:ntyp1) + common/fourier/ b1(2,maxres),b2(2,maxres),b(13,-ntyp:ntyp), + & bnew1(3,2,-ntyp:ntyp),bnew2(3,2,-ntyp:ntyp), + & cc(2,2,-ntyp:ntyp), + & dd(2,2,-ntyp:ntyp),eeold(2,2,-ntyp:ntyp), + & eenew(2,-ntyp:ntyp), & ee(2,2,maxres), - & ctilde(2,2,-maxtor:maxtor), - & dtilde(2,2,-maxtor:maxtor),b1tilde(2,0:maxres), + & ctilde(2,2,-ntyp:ntyp), + & dtilde(2,2,-ntyp:ntyp),b1tilde(2,maxres), & b2tilde(2,maxres), & gtb1(2,maxres),gtb2(2,maxres),gtEE(2,2,maxres), - & nloctyp - + & nloctyp,iloctyp,itype2loc diff --git a/source/unres/src_MD-M/Makefile b/source/unres/src_MD-M/Makefile deleted file mode 100644 index f83b04c..0000000 --- a/source/unres/src_MD-M/Makefile +++ /dev/null @@ -1,130 +0,0 @@ -#source /opt/pgi/linux86-64/13.7/mpi.csh -################################################################### - - -FC= mpif90 -OPT = -Minfo - -FFLAGS = -c ${OPT} -FFLAGS1 = -c -g -FFLAGS2 = -c -g -O0 -FFLAGSE = -c -fast -Minline=name:scalar2,scalar,transpose2,matvec2,prodmat3 -Minfo - - -LIBS = xdrf/libxdrf.a - -ARCH = LINUX -PP = /lib/cpp -P - - -all: no_option - @echo "give optin GAB or E0LL2Y" - -.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 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 ssMD.o isnan.o permut.o together.o - -no_option: - -GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI \ - -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -GAB: BIN = ../../../bin/unres/MD/unres_pgf90_mpi_GAB.exe -GAB: ${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 -DUNRES -DISNAN -DMP -DMPI \ - -DSPLITELE -DLANG0 -E0LL2Y: BIN = ../../../bin/unres/MD/unres_pgf90_mpi_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 - -isnan.o: isnan.f - ${FC} -Kieee -c isnan.f - 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/Makefile_MPICH_ifort b/source/unres/src_MD-M/Makefile_MPICH_ifort index 0d2a6d4..7ce0741 100644 --- a/source/unres/src_MD-M/Makefile_MPICH_ifort +++ b/source/unres/src_MD-M/Makefile_MPICH_ifort @@ -5,11 +5,13 @@ INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh FC= ifort OPT = -O3 -ip +OPT = -g -CA -CB FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include FFLAGS1 = -c -g -CA -CB -I$(INSTALL_DIR)/include FFLAGS2 = -c -g -O0 -I$(INSTALL_DIR)/include FFLAGSE = -c -O3 -ipo -opt_report -I$(INSTALL_DIR)/include +FFLAGSE = ${FFLAGS} LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a @@ -47,7 +49,7 @@ no_option: GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -GAB: BIN = ../../../bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe +GAB: BIN = /users/adam/bin/unres/MD/unres-mult-symetr_KCC_ifort_MPICH_GAB.exe GAB: ${object} xdrf/libxdrf.a cc -o compinfo compinfo.c ./compinfo | true @@ -56,7 +58,7 @@ GAB: ${object} xdrf/libxdrf.a 4P: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -4P: BIN = ../../../bin/unres/MD/unres-mult-symetr_ifort_MPICH_4P.exe +4P: BIN = /users/adam/bin/unres-mult-symetr_KCC_ifort_MPICH_4P.exe 4P: ${object} xdrf/libxdrf.a cc -o compinfo compinfo.c ./compinfo | true @@ -65,13 +67,22 @@ GAB: ${object} xdrf/libxdrf.a E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ -DSPLITELE -DLANG0 -E0LL2Y: BIN = ../../../bin/unres/MD/unres-mult-symetr_ifort_MPICH_E0LL2Y.exe +E0LL2Y: BIN = /users/adam/bin/unres-mult-symetr_KCC_ifort_MPICH_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} +NEWCORR: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ + -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR: BIN = /users/adam/bin/unres-mult-symetr_KCC_ifort_MPICH_E0LL2Y-NEWCORR.exe +NEWCORR: ${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 diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F index 270f4cc..80dd213 100644 --- a/source/unres/src_MD-M/checkder_p.F +++ b/source/unres/src_MD-M/checkder_p.F @@ -2,6 +2,7 @@ C Check the gradient of Cartesian coordinates in internal coordinates. implicit real*8 (a-h,o-z) include 'DIMENSIONS' + include 'COMMON.CONTROL' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CHAIN' @@ -14,8 +15,9 @@ C Check the gradient of Cartesian coordinates in internal coordinates. * Check the gradient of the virtual-bond and SC vectors in the internal * coordinates. * - aincr=1.5d-7 - aincr2=2.5d-8 + print '("Calling CHECK_ECART",1pd12.3)',aincr + write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr + aincr2=0.5d0*aincr call cartder write (iout,'(a)') '**************** dx/dalpha' write (iout,'(a)') @@ -175,6 +177,7 @@ C---------------------------------------------------------------------------- C Check the gradient of the energy in Cartesian coordinates. implicit real*8 (a-h,o-z) include 'DIMENSIONS' + include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.IOUNITS' @@ -191,8 +194,8 @@ C Check the gradient of the energy in Cartesian coordinates. nf=0 nfl=0 call zerograd - aincr=1.0D-7 - print '(a)','CG processor',me,' calling CHECK_CART.' + print '("Calling CHECK_ECART",1pd12.3)',aincr + write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr nf=0 icall=0 call geom_to_var(nvar,x) @@ -282,8 +285,10 @@ c rlambd=0.3d0 c call intcartderiv c call checkintcartgrad call zerograd - aincr=8.0D-7 - write(iout,*) 'Calling CHECK_ECARTINT.' +c aincr=8.0D-7 +c aincr=1.0D-7 + print '("Calling CHECK_ECARTINT",1pd12.3)',aincr + write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr nf=0 icall=0 call geom_to_var(nvar,x) @@ -298,8 +303,10 @@ c call checkintcartgrad write (iout,*) "exit cartgrad" call flush(iout) icall =1 + write (iout,*) "gcard and gxcart" do i=1,nres - write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) + write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), + & (gxcart(j,i),j=1,3) enddo do j=1,3 grad_s(j,0)=gcart(j,0) @@ -628,6 +635,7 @@ c---------------------------------------------------------------------------- C Check the gradient of energy in internal coordinates. implicit real*8 (a-h,o-z) include 'DIMENSIONS' + include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.IOUNITS' @@ -642,8 +650,9 @@ C Check the gradient of energy in internal coordinates. character*6 key external fdum call zerograd - aincr=1.0D-7 - print '(a)','Calling CHECK_INT.' +c aincr=1.0D-7 + print '("Calling CHECK_INT",1pd12.3)',aincr + write (iout,'("Calling CHECK_INT",1pd12.3)') aincr nf=0 nfl=0 icg=1 @@ -674,12 +683,12 @@ cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar) xi=x(i) x(i)=xi-0.5D0*aincr call var_to_geom(nvar,x) - call chainbuild + call chainbuild_extconf call etotal(energia1(0)) etot1=energia1(0) x(i)=xi+0.5D0*aincr call var_to_geom(nvar,x) - call chainbuild + call chainbuild_extconf call etotal(energia2(0)) etot2=energia2(0) gg(i)=(etot2-etot1)/aincr diff --git a/source/unres/src_MD-M/cinfo.f b/source/unres/src_MD-M/cinfo.f index 9d0e3bd..d77b707 100644 --- a/source/unres/src_MD-M/cinfo.f +++ b/source/unres/src_MD-M/cinfo.f @@ -1,30 +1,37 @@ C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C -C 0 40360 18 +C 0 40360 71 subroutine cinfo include 'COMMON.IOUNITS' write(iout,*)'++++ Compile info ++++' - write(iout,*)'Version 0.40360 build 18' - write(iout,*)'compiled Fri Nov 27 10:42:05 2015' - write(iout,*)'compiled by adasko@cuda1' + write(iout,*)'Version 0.40360 build 71' + write(iout,*)'compiled Sun Jan 3 21:03:16 2016' + write(iout,*)'compiled by adam@piasek4' write(iout,*)'OS name: Linux ' - write(iout,*)'OS release: 3.2.0-32-generic ' + write(iout,*)'OS release: 3.2.0-91-generic ' write(iout,*)'OS version:', - & ' #51-Ubuntu SMP Wed Sep 26 21:33:09 UTC 2012 ' + & ' #129-Ubuntu SMP Wed Sep 9 10:56:06 UTC 2015 ' write(iout,*)'flags:' - write(iout,*)'FC= mpif90' - write(iout,*)'OPT = -Minfo' - write(iout,*)'FFLAGS = -c ${OPT} ' - write(iout,*)'FFLAGS1 = -c -g ' - write(iout,*)'FFLAGS2 = -c -g -O0 ' - write(iout,*)'FFLAGSE = -c -fast -Minline=name:scalar2,scala...' - write(iout,*)'LIBS = xdrf/libxdrf.a' + write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...' + write(iout,*)'FC= ifort' + write(iout,*)'OPT = -O3 -ip ' + write(iout,*)'OPT = -g -CA -CB' + write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include ' + write(iout,*)'FFLAGS1 = -c -g -CA -CB -I$(INSTALL_DIR)/inclu...' + write(iout,*)'FFLAGS2 = -c -g -O0 -I$(INSTALL_DIR)/include ' + write(iout,*)'FFLAGSE = -c -O3 -ipo -opt_report -I$(INSTALL...' + write(iout,*)'FFLAGSE = ${FFLAGS}' + write(iout,*)'LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdr...' write(iout,*)'ARCH = LINUX' write(iout,*)'PP = /lib/cpp -P' write(iout,*)'object = unres.o arcos.o cartprint.o chainbuild...' - write(iout,*)'GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES ...' - write(iout,*)'GAB: BIN = ../../../bin/unres/MD/unres_pgf90_mp...' - write(iout,*)'E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNR...' - write(iout,*)'E0LL2Y: BIN = ../../../bin/unres/MD/unres_pgf90...' + write(iout,*)'GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 ...' + write(iout,*)'GAB: BIN = /users/adam/bin/unres/MD/unres-mult-...' + write(iout,*)'4P: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -D...' + write(iout,*)'4P: BIN = /users/adam/bin/unres-mult-symetr_KCC...' + write(iout,*)'E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD...' + write(iout,*)'E0LL2Y: BIN = /users/adam/bin/unres-mult-symetr...' + write(iout,*)'NEWCORR: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAM...' + write(iout,*)'NEWCORR: BIN = /users/adam/bin/unres-mult-symet...' write(iout,*)'++++ End of compile info ++++' 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 f524af3..8a2d035 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2783,28 +2783,28 @@ c write(iout,*) 'nphi=',nphi,nres #endif #ifdef NEWCORR if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itype2loc(itype(i-2)) else - iti=ntortyp+1 + iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp+1 + iti1=nloctyp endif c write(iout,*),i - b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0) + b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0d0) & +bnew1(2,1,iti)*dsin(theta(i-1)) - & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0) + & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0d0) gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 & +bnew1(2,1,iti)*dcos(theta(i-1)) & -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) - b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0) + b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0d0) & +bnew2(2,1,iti)*dsin(theta(i-1)) - & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0) + & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0d0) c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 @@ -2843,15 +2843,15 @@ c write (iout,*) 'theta=', theta(i-1) enddo #else if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itype2loc(itype(i-2)) else - iti=ntortyp+1 + iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp+1 + iti1=nloctyp endif b1(1,i-2)=b(3,iti) b1(2,i-2)=b(5,iti) @@ -2940,15 +2940,15 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then endif c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itype2loc(itype(i-2)) else - iti=ntortyp + iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp + iti1=nloctyp endif cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) @@ -2961,8 +2961,8 @@ c if (i .gt. iatel_s+2) then call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2)) c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" #endif -c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti), -c & EE(1,2,iti),EE(2,2,iti) +c write(iout,*) "co jest kurwa", iti, EE(1,1,i),EE(2,1,i), +c & EE(1,2,iti),EE(2,2,i) call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2)) call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2)) c write(iout,*) "Macierz EUG", @@ -2997,18 +2997,24 @@ c & eug(2,2,i-2) c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then if (itype(i-1).le.ntyp) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp + iti1=nloctyp endif else - iti1=ntortyp + iti1=nloctyp endif do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo -C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) -c write (iout,*) 'mu ',mu(:,i-2),i-2 +#ifdef MUOUT + write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1), + & rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2), + & -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2), + & dsqrt(b2(1,i-1)**2+b2(2,i-1)**2) + & +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2), + & ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itype2loc(iti)) +#endif cd write (iout,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) @@ -3295,11 +3301,11 @@ c endif #endif #endif cd do i=1,nres -cd iti = itortyp(itype(i)) +cd iti = itype2loc(itype(i)) cd write (iout,*) i cd do j=1,2 cd write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)') -cd & (EE(j,k,iti),k=1,2),(Ug(j,k,i),k=1,2),(EUg(j,k,i),k=1,2) +cd & (EE(j,k,i),k=1,2),(Ug(j,k,i),k=1,2),(EUg(j,k,i),k=1,2) cd enddo cd enddo return @@ -3417,21 +3423,23 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups C C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end - if (i.le.1) cycle +c if (i.le.1) cycle C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+4).gt.nres) - & .or.((i-1).le.0) +C Adam: Unnecessary: handled by iturn3_end and iturn3_start +c & .or.((i+4).gt.nres) +c & .or.((i-1).le.0) C end of changes by Ana & .or. itype(i+2).eq.ntyp1 & .or. itype(i+3).eq.ntyp1) cycle - if(i.gt.1)then - if(itype(i-1).eq.ntyp1)cycle - end if - if(i.LT.nres-3)then - if (itype(i+4).eq.ntyp1) cycle - end if +C Adam: Instructions below will switch off existing interactions +c if(i.gt.1)then +c if(itype(i-1).eq.ntyp1)cycle +c end if +c if(i.LT.nres-3)then +c if (itype(i+4).eq.ntyp1) cycle +c end if dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3456,14 +3464,14 @@ C end of changes by Ana if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+5).gt.nres) - & .or.((i-1).le.0) +c & .or.((i+5).gt.nres) +c & .or.((i-1).le.0) C end of changes suggested by Ana & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 - & .or. itype(i).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +c & .or. itype(i+5).eq.ntyp1 +c & .or. itype(i).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3523,14 +3531,14 @@ c CTU KURWA do i=iatel_s,iatel_e C do i=75,75 - if (i.le.1) cycle +c if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+2).gt.nres) - & .or.((i-1).le.0) +c & .or.((i+2).gt.nres) +c & .or.((i-1).le.0) C end of changes by Ana - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +c & .or. itype(i+2).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3586,11 +3594,11 @@ C write (iout,*) i,j if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((j+2).gt.nres) - & .or.((j-1).le.0) +c & .or.((j+2).gt.nres) +c & .or.((j-1).le.0) C end of changes by Ana - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +c & .or.itype(j+2).eq.ntyp1 +c & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j @@ -4832,9 +4840,9 @@ c write(iout,*)"WCHODZE W PROGRAM" a_temp(1,2)=a23 a_temp(2,1)=a32 a_temp(2,2)=a33 - iti1=itortyp(itype(i+1)) - iti2=itortyp(itype(i+2)) - iti3=itortyp(itype(i+3)) + iti1=itype2loc(itype(i+1)) + iti2=itype2loc(itype(i+2)) + iti3=itype2loc(itype(i+3)) c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) @@ -7473,11 +7481,12 @@ C The rigorous attempt to derive energy function include 'COMMON.TORCNSTR' include 'COMMON.CONTROL' logical lprn - double precision thybt1(maxtermkcc),thybt2(maxtermkcc) +c double precision thybt1(maxtermkcc),thybt2(maxtermkcc) C Set lprn=.true. for debugging lprn=.false. c lprn=.true. C print *,"wchodze kcc" + if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode if (tor_mode.ne.2) then etors=0.0D0 endif @@ -7497,60 +7506,86 @@ c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle sumnonchebyshev=0.0d0 sumchebyshev=0.0d0 C to avoid multiple devision by 2 - theti22=0.5d0*theta(i) +c theti22=0.5d0*theta(i) C theta 12 is the theta_1 /2 C theta 22 is theta_2 /2 - theti12=0.5d0*theta(i-1) +c theti12=0.5d0*theta(i-1) C and appropriate sinus function - sinthet2=dsin(theta(i)) sinthet1=dsin(theta(i-1)) + sinthet2=dsin(theta(i)) costhet1=dcos(theta(i-1)) costhet2=dcos(theta(i)) +c Cosines of halves thetas + costheti12=0.5d0*(1.0d0+costhet1) + costheti22=0.5d0*(1.0d0+costhet2) C to speed up lets store its mutliplication - sint1t2=sinthet2*sinthet1 + sint1t2=sinthet2*sinthet1 + sint1t2n=1.0d0 C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma) C +d_n*sin(n*gamma)) * C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2))) C we have two sum 1) Non-Chebyshev which is with n and gamma + etori=0.0d0 do j=1,nterm_kcc(itori,itori1) + nval=nterm_kcc_Tb(itori,itori1) v1ij=v1_kcc(j,itori,itori1) v2ij=v2_kcc(j,itori,itori1) +c write (iout,*) "i",i," j",j," v1",v1ij," v2",v2ij C v1ij is c_n and d_n in euation above cosphi=dcos(j*phii) sinphi=dsin(j*phii) - sint1t2n=sint1t2**j - sumnonchebyshev= - & sint1t2n*(v1ij*cosphi+v2ij*sinphi) - actval=sint1t2n*(v1ij*cosphi+v2ij*sinphi) + sint1t2n1=sint1t2n + sint1t2n=sint1t2n*sint1t2 + sumth1tyb1=tschebyshev(1,nval,v11_chyb(1,j,itori,itori1), + & costheti12) + gradth1tyb1=-0.5d0*sinthet1*gradtschebyshev(0,nval-1, + & v11_chyb(1,j,itori,itori1),costheti12) +c write (iout,*) "v11",(v11_chyb(k,j,itori,itori1),k=1,nval), +c & " sumth1tyb1",sumth1tyb1," gradth1tyb1",gradth1tyb1 + sumth2tyb1=tschebyshev(1,nval,v21_chyb(1,j,itori,itori1), + & costheti22) + gradth2tyb1=-0.5d0*sinthet2*gradtschebyshev(0,nval-1, + & v21_chyb(1,j,itori,itori1),costheti22) +c write (iout,*) "v21",(v21_chyb(k,j,itori,itori1),k=1,nval), +c & " sumth2tyb1",sumth2tyb1," gradth2tyb1",gradth2tyb1 + sumth1tyb2=tschebyshev(1,nval,v12_chyb(1,j,itori,itori1), + & costheti12) + gradth1tyb2=-0.5d0*sinthet1*gradtschebyshev(0,nval-1, + & v12_chyb(1,j,itori,itori1),costheti12) +c write (iout,*) "v12",(v12_chyb(k,j,itori,itori1),k=1,nval), +c & " sumth1tyb2",sumth1tyb2," gradth1tyb2",gradth1tyb2 + sumth2tyb2=tschebyshev(1,nval,v22_chyb(1,j,itori,itori1), + & costheti22) + gradth2tyb2=-0.5d0*sinthet2*gradtschebyshev(0,nval-1, + & v22_chyb(1,j,itori,itori1),costheti22) +c write (iout,*) "v22",(v22_chyb(k,j,itori,itori1),k=1,nval), +c & " sumth2tyb2",sumth2tyb2," gradth2tyb2",gradth2tyb2 C etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi) C if (energy_dec) etors_ii=etors_ii+ C & v1ij*cosphi+v2ij*sinphi C glocig is the gradient local i site in gamma - glocig=j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n + actval1=v1ij*cosphi*(1.0d0+sumth1tyb1+sumth2tyb1) + actval2=v2ij*sinphi*(1.0d0+sumth1tyb2+sumth2tyb2) + etori=etori+sint1t2n*(actval1+actval2) + glocig=glocig+ + & j*sint1t2n*(v2ij*cosphi*(1.0d0+sumth1tyb2+sumth2tyb2) + & -v1ij*sinphi*(1.0d0+sumth1tyb1+sumth2tyb1)) C now gradient over theta_1 - glocit1=actval/sinthet1*j*costhet1 - glocit2=actval/sinthet2*j*costhet2 + glocit1=glocit1+ + & j*sint1t2n1*costhet1*sinthet2*(actval1+actval2)+ + & sint1t2n*(v1ij*cosphi*gradth1tyb1+v2ij*sinphi*gradth1tyb2) + glocit2=glocit2+ + & j*sint1t2n1*sinthet1*costhet2*(actval1+actval2)+ + & sint1t2n*(v1ij*cosphi*gradth2tyb1+v2ij*sinphi*gradth2tyb2) C now the Czebyshev polinominal sum - do k=1,nterm_kcc_Tb(itori,itori1) - thybt1(k)=v1_chyb(k,j,itori,itori1) - thybt2(k)=v2_chyb(k,j,itori,itori1) +c do k=1,nterm_kcc_Tb(itori,itori1) +c thybt1(k)=v1_chyb(k,j,itori,itori1) +c thybt2(k)=v2_chyb(k,j,itori,itori1) C thybt1(k)=0.0 C thybt2(k)=0.0 - enddo - sumth1thyb=tschebyshev - & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theti12)**2) - gradthybt1=gradtschebyshev - & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1), - & dcos(theti12)**2) - & *dcos(theti12)*(-dsin(theti12)) - sumth2thyb=tschebyshev - & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theti22)**2) - gradthybt2=gradtschebyshev - & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1), - & dcos(theti22)**2) - & *dcos(theti22)*(-dsin(theti22)) +c enddo C print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2, C & gradtschebyshev C & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1), @@ -7558,22 +7593,19 @@ C & dcos(theti22)**2), C & dsin(theti22) C now overal sumation - etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev C print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb + enddo ! j + etors=etors+etori C derivative over gamma - gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig - & *(1.0d0+sumth1thyb+sumth2thyb) + gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig C derivative over theta1 - gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor* - & (glocit1*(1.0d0+sumth1thyb+sumth2thyb)+ - & sumnonchebyshev*gradthybt1) + gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*glocit1 C now derivative over theta2 - gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor* - & (glocit2*(1.0d0+sumth1thyb+sumth2thyb)+ - & sumnonchebyshev*gradthybt2) - enddo + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2 + if (lprn) + & write (iout,*) i-2,i-1,itype(i-2),itype(i-1),itori,itori1, + & theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori enddo - C gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci ! 6/20/98 - dihedral angle constraints if (tor_mode.ne.2) then @@ -7622,7 +7654,8 @@ C Set lprn=.true. for debugging lprn=.false. c lprn=.true. C print *,"wchodze kcc" - if (tormode.ne.2) etheta=0.0D0 + if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode + if (tor_mode.ne.2) etheta=0.0D0 do i=ithet_start,ithet_end c print *,i,itype(i-1),itype(i),itype(i-2) if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 @@ -7635,6 +7668,8 @@ c print *,i,itype(i-1),itype(i),itype(i-2) enddo sumth1thyb=tschebyshev & (1,nbend_kcc_Tb(iti),thybt1(1),costhet) + if (lprn) write (iout,*) i-1,itype(i-1),iti,theta(i)*rad2deg, + & sumth1thyb ihelp=nbend_kcc_Tb(iti)-1 gradthybt1=gradtschebyshev & (0,ihelp,thybt1(1),costhet) @@ -7643,7 +7678,7 @@ C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0) gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang* & gradthybt1*sinthet*(-0.5d0) enddo - if (tormode.ne.2) then + if (tor_mode.ne.2) then ethetacnstr=0.0d0 C print *,ithetaconstr_start,ithetaconstr_end,"TU" do i=ithetaconstr_start,ithetaconstr_end @@ -8839,9 +8874,9 @@ C--------------------------------------------------------------------------- & auxmat(2,2) iti1 = itortyp(itype(i+1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1)) + itj1 = itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -8929,16 +8964,16 @@ cd write (iout,*) "a_chujkl",((a_chuj(iii,jjj,kk,k),iii=1,2),jjj=1,2) if (l.eq.j+1) then C parallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itype2loc(itype(i)) else - iti=ntortyp + iti=nloctyp endif - itk1=itortyp(itype(k+1)) - itj=itortyp(itype(j)) + itk1=itype2loc(itype(k+1)) + itj=itype2loc(itype(j)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif C A1 kernel(j+1) A2T cd do iii=1,2 @@ -9082,17 +9117,17 @@ C End vectors else C Antiparallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itype2loc(itype(i)) else - iti=ntortyp + iti=nloctyp endif - itk1=itortyp(itype(k+1)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + itk1=itype2loc(itype(k+1)) + itl=itype2loc(itype(l)) + itj=itype2loc(itype(j)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif C A2 kernel(j-1)T A1T call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), @@ -9430,9 +9465,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)) + itk=itype2loc(itype(k)) + itl=itype2loc(itype(l)) + itj=itype2loc(itype(j)) eello5_1=0.0d0 eello5_2=0.0d0 eello5_3=0.0d0 @@ -9501,7 +9536,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) @@ -9582,7 +9617,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) @@ -9655,7 +9690,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) @@ -9952,7 +9987,7 @@ C o o o o C C i i C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - itk=itortyp(itype(k)) + itk=itype2loc(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)) @@ -10244,16 +10279,16 @@ 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)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) + itk=itype2loc(itype(k)) + itk1=itype2loc(itype(k+1)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) @@ -10262,7 +10297,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) @@ -10360,24 +10395,24 @@ 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)) + iti=itype2loc(itype(i)) + itj=itype2loc(itype(j)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif - itk=itortyp(itype(k)) + itk=itype2loc(itype(k)) if (k.lt.nres-1) then - itk1=itortyp(itype(k+1)) + itk1=itype2loc(itype(k+1)) else - itk1=ntortyp + itk1=nloctyp endif - itl=itortyp(itype(l)) + itl=itype2loc(itype(l)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp 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, @@ -10597,11 +10632,11 @@ c j=i+4 k=i+1 l=i+3 - iti=itortyp(itype(i)) - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + iti=itype2loc(itype(i)) + itk=itype2loc(itype(k)) + itk1=itype2loc(itype(k+1)) + itl=itype2loc(itype(l)) + itj=itype2loc(itype(j)) cd write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj cd write (2,*) 'i',i,' k',k,' j',j,' l',l cd if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then @@ -11441,7 +11476,7 @@ C-------------------------------------------------------------------------- include "DIMENSIONS" integer i,m,n double precision x(n+1),y,yy(0:maxvar),aux -c Tschebyshev polynomial. Note that the first term is omitted +c Tschebyshev polynomial. Note that the first term is omitted c m=0: the constant term is included c m=1: the constant term is not included yy(0)=1.0d0 diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F index 4861790..4cfa18b 100644 --- a/source/unres/src_MD-M/gradient_p.F +++ b/source/unres/src_MD-M/gradient_p.F @@ -347,6 +347,7 @@ C------------------------------------------------------------------------- include 'COMMON.MD' include 'COMMON.SCCOR' include 'COMMON.SHIELD' + maxshieldlist=0 C C Initialize Cartesian-coordinate gradient C diff --git a/source/unres/src_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F index 562ea78..227a119 100644 --- a/source/unres/src_MD-M/intcartderiv.F +++ b/source/unres/src_MD-M/intcartderiv.F @@ -119,7 +119,7 @@ c the conventional case c Obtaining the gamma derivatives from sine derivative if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. & phi(i).gt.pi34.and.phi(i).le.pi.or. - & phi(i).gt.-pi.and.phi(i).le.-pi34) then + & phi(i).ge.-pi.and.phi(i).le.-pi34) then call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1) call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2) call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3) @@ -199,7 +199,7 @@ cc write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4 c Obtaining the gamma derivatives from sine derivative if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or. & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or. - & tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then + & tauangle(1,i).ge.-pi.and.tauangle(1,i).le.-pi34) then call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1) call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2) call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3) @@ -346,7 +346,7 @@ c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) c Obtaining the gamma derivatives from sine derivative if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or. & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or. - & tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then + & tauangle(3,i).ge.-pi.and.tauangle(3,i).le.-pi34) then call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1) call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2) call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3) @@ -425,7 +425,7 @@ c write (iout,*) "i",i," cosa",cosa," sina",sina," sino",sino c obtaining the derivatives of omega from sines if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or. & omeg(i).gt.pi34.and.omeg(i).le.pi.or. - & omeg(i).gt.-pi.and.omeg(i).le.-pi34) then + & omeg(i).ge.-pi.and.omeg(i).le.-pi34) then fac15=dcos(theta(i+1))/(dsin(theta(i+1))* & dsin(theta(i+1))) fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i))) diff --git a/source/unres/src_MD-M/minimize_p.F b/source/unres/src_MD-M/minimize_p.F index 326c02d..ef24c0c 100644 --- a/source/unres/src_MD-M/minimize_p.F +++ b/source/unres/src_MD-M/minimize_p.F @@ -263,7 +263,7 @@ c endif cd print *,'func',nf,nfl,icg call var_to_geom(n,x) call zerograd - call chainbuild + call chainbuild_extconf cd write (iout,*) 'ETOTAL called from FUNC' call etotal(energia(0)) call sum_gradient @@ -297,7 +297,7 @@ c endif icg=mod(nf,2)+1 call var_to_geom_restr(n,x) call zerograd - call chainbuild + call chainbuild_extconf cd write (iout,*) 'ETOTAL called from FUNC' call etotal(energia(0)) call sum_gradient diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 438f544..e96dfc3 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -101,6 +101,8 @@ c enddo endif C reading lipid parameters + write (iout,*) "iliptranpar",iliptranpar + call flush(iout) read(iliptranpar,*) pepliptran do i=1,ntyp read(iliptranpar,*) liptranene(i) @@ -219,6 +221,8 @@ C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 C read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2, & ntheterm3,nsingle,ndouble + write (iout,*) "ithep",ithep + call flush(iout) nntheterm=max0(ntheterm,ntheterm2,ntheterm3) read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1) do i=-ntyp1,-1 @@ -361,11 +365,12 @@ C Control printout of the coefficients of virtual-bond-angle potentials C if (lprint) then write (iout,'(//a)') 'Parameter of virtual-bond-angle potential' - do i=1,nthetyp+1 - do j=1,nthetyp+1 - do k=1,nthetyp+1 + do iblock=1,2 + do i=0,nthetyp + do j=-nthetyp,nthetyp + do k=-nthetyp,nthetyp write (iout,'(//4a)') - & 'Type ',onelett(i),onelett(j),onelett(k) + & 'Type ',toronelet(i),toronelet(j),toronelet(k) write (iout,'(//a,10x,a)') " l","a[l]" write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock) write (iout,'(i2,1pe15.5)') @@ -394,6 +399,7 @@ C enddo enddo enddo + enddo enddo call flush(iout) endif @@ -661,8 +667,9 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) close (itorp) if (lprint) then write (iout,'(/a/)') 'Torsional constants:' - do i=1,ntortyp - do j=1,ntortyp + do iblock=1,2 + do i=0,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 write (iout,*) 'ityp',i,' jtyp',j write (iout,*) 'Fourier constants' do k=1,nterm(i,j,iblock) @@ -676,6 +683,7 @@ c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo enddo enddo + enddo endif C @@ -776,42 +784,79 @@ C Martix of D parameters for two dimesional fourier series endif #endif C read Czybyshev torsional parameters - read (itorkcc,*,end=121,err=121) nkcctyp - read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp) - do i=-ntyp,-1 + read (itorkcc,*,end=121,err=121) nkcctyp + read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp) + do i=-ntyp,-1 itortyp_kcc(i)=-itortyp_kcc(-i) - enddo - do i=1,nkcctyp - do j=1,nkcctyp + enddo + do i=0,nkcctyp + do j=0,nkcctyp C first we read the cos and sin gamma parameters read (itorkcc,*,end=121,err=121) & nterm_kcc(j,i),nterm_kcc_Tb(j,i) C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) - do k=1,nterm_kcc(j,i) - do l=1,nterm_kcc_Tb(j,i) - read (itorkcc,*,end=121,err=121) - & v1_chyb(l,k,j,i) - enddo - do l=1,nterm_kcc_Tb(j,i) - read (itorkcc,*,end=121,err=121) - & v2_chyb(l,k,j,i) - enddo - read (itorkcc,*,end=121,err=121) - & v1_kcc(k,j,i) - read (itorkcc,*,end=121,err=121) - & v2_kcc(k,j,i) - enddo + do k=1,nterm_kcc(j,i) + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) v11_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i) + enddo + do l=1,nterm_kcc_Tb(j,i) + read (itorkcc,*,end=121,err=121) v22_chyb(l,k,j,i) + enddo + read (itorkcc,*,end=121,err=121) v1_kcc(k,j,i) + read (itorkcc,*,end=121,err=121) v2_kcc(k,j,i) + enddo + enddo + enddo + if (lprint) then +c Print valence-torsional parameters + write (iout,'(a)') + & "Parameters of the valence-torsional potentials" + do i=0,nkcctyp + do j=0,nkcctyp + write (iout,'(3a)') "Type ",toronelet(i),toronelet(j) + write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb" + do k=1,nterm_kcc(j,i) + write (iout,'(i5,f15.10,i5,2f15.10)') + & k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i) + do l=2,nterm_kcc_Tb(j,i) + write (iout,'(20x,i5,2f15.10)') + & l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i) enddo - enddo + write (iout,'(i5,f15.10,i5,2f15.10)') + & k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i) + do l=2,nterm_kcc_Tb(j,i) + write (iout,'(20x,i5,2f15.10)') + & l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i) + enddo + write (iout,'(a)') + enddo + enddo + enddo + endif C here will be the apropriate recalibrating for D-aminoacid C read (ithetkcc,*,end=121,err=121) nkcctyp - do i=1,nkcctyp + do i=0,nkcctyp read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i) - do j=1,nbend_kcc_Tb(i) - read (ithetkcc,*,end=121,err=121) - & v1bend_chyb(j,i) - enddo + do j=1,nbend_kcc_Tb(i) + read (ithetkcc,*,end=121,err=121) v1bend_chyb(j,i) + enddo + enddo + if (lprint) then + write (iout,'(a)') + & "Parameters of the valence-only potentials" + do i=0,nkcctyp + write (iout,'(2a)') "Type ",toronelet(i) + do k=1,nbend_kcc_Tb(i) + write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i) enddo + enddo + endif C Read of Side-chain backbone correlation parameters C Modified 11 May 2012 by Adasko CCC @@ -919,7 +964,7 @@ 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 @@ -928,6 +973,7 @@ cc maxinter is maximum interaction sites #endif if (lprint) then write (iout,'(/a/)') 'Torsional constants:' + do l=1,maxinter do i=1,nsccortyp do j=1,nsccortyp write (iout,*) 'ityp',i,' jtyp',j @@ -942,6 +988,7 @@ cc maxinter is maximum interaction sites enddo enddo enddo + enddo endif C @@ -953,7 +1000,29 @@ C write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp - +#ifdef NEWCORR + read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp) + read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1) + itype2loc(ntyp1)=nloctyp + iloctyp(nloctyp)=ntyp1 +#else + do i=1,ntyp1 + itype2loc(i)=itortyp(i) + enddo + iloctyp(0)=10 + iloctyp(1)=9 + iloctyp(2)=20 + iloctyp(3)=ntyp1 +#endif + do i=1,ntyp1 + itype2loc(-i)=-itype2loc(i) + enddo + do i=1,nloctyp + iloctyp(-i)=-iloctyp(i) + enddo + write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) + write (iout,*) "nloctyp",nloctyp, + & " iloctyp",(iloctyp(i),i=0,nloctyp) do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) @@ -984,14 +1053,14 @@ c B2(1,i) = b(2) c B2(2,i) = b(4) c B2(1,-i) =b(2) c B2(2,-i) =-b(4) - B1tilde(1,i) = b(3,i) - B1tilde(2,i) =-b(5,i) +cc B1tilde(1,i) = b(3,i) +cc B1tilde(2,i) =-b(5,i) C B1tilde(1,-i) =-b(3,i) C B1tilde(2,-i) =b(5,i) - b1tilde(1,i)=0.0d0 - b1tilde(2,i)=0.0d0 - B2(1,i) = b(2,i) - B2(2,i) = b(4,i) +cc b1tilde(1,i)=0.0d0 +cc b1tilde(2,i)=0.0d0 +cc B2(1,i) = b(2,i) +cc B2(2,i) = b(4,i) C B2(1,-i) =b(2,i) C B2(2,-i) =-b(4,i) diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 9f72992..c691135 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -179,6 +179,7 @@ C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then else icheckgrad=3 endif + call reada(controlcard,'DELTA',aincr,1.0d-7) elseif (index(controlcard,'THREAD').gt.0) then modecalc=2 call readi(controlcard,'THREAD',nthread,0)