From f9783c27de185ff3e1d6e874a37a1a8c4d3c93a3 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Sun, 7 Jun 2020 12:15:53 +0200 Subject: [PATCH] Adam's lipid and dfa corrections --- source/cluster/wham/src-HCD/Makefile | 2 +- source/cluster/wham/src-HCD/Makefile-MPICH-ifort | 84 ++- .../wham/src-HCD/Makefile-MPICH-ifort-prometheus | 82 ++- source/cluster/wham/src-HCD/Makefile-tryton | 2 +- source/cluster/wham/src-HCD/energy_p_new.F | 128 ++-- .../wham/src-HCD/include_unres/COMMON.INTERACT | 7 +- source/cluster/wham/src-HCD/initialize.f | 99 --- source/cluster/wham/src-HCD/readrtns.F | 1 + source/cluster/wham/src-HCD/sizesclu.dat | 2 +- source/unres/src-HCD-5D/COMMON.DFA | 21 +- source/unres/src-HCD-5D/COMMON.INTERACT | 4 +- source/unres/src-HCD-5D/DIMENSIONS | 4 +- source/unres/src-HCD-5D/MD_A-MTS.F | 1 + source/unres/src-HCD-5D/checkder_p.F | 9 + source/unres/src-HCD-5D/dfa.F | 52 +- source/unres/src-HCD-5D/energy_p_new-sep_barrier.F | 322 +++++++--- source/unres/src-HCD-5D/energy_p_new_barrier.F | 308 +++++++--- source/unres/src-HCD-5D/energy_split-sep.F | 8 +- source/unres/src-HCD-5D/gradient_p.F | 9 +- source/unres/src-HCD-5D/parmread.F | 13 +- source/unres/src-HCD-5D/readpdb-mult.F | 1 + source/unres/src-HCD-5D/readpdb.F | 631 -------------------- source/unres/src-HCD-5D/readrtns_CSA.F | 7 +- source/unres/src-HCD-5D/unres.F | 34 +- source/wham/src-HCD/Makefile | 2 +- source/wham/src-HCD/Makefile_MPICH_ifort | 100 +++- source/wham/src-HCD/energy_p_new.F | 152 +++-- source/wham/src-HCD/include_unres/COMMON.INTERACT | 6 +- source/wham/src-HCD/initialize_p.F | 2 +- source/wham/src-HCD/oligomer.F | 4 +- source/wham/src-HCD/parmread.F | 2 + source/wham/src-HCD/proc_cont.f | 4 +- 32 files changed, 1008 insertions(+), 1095 deletions(-) delete mode 100644 source/cluster/wham/src-HCD/initialize.f delete mode 100644 source/unres/src-HCD-5D/readpdb.F diff --git a/source/cluster/wham/src-HCD/Makefile b/source/cluster/wham/src-HCD/Makefile index 8aee570..693492e 120000 --- a/source/cluster/wham/src-HCD/Makefile +++ b/source/cluster/wham/src-HCD/Makefile @@ -1 +1 @@ -Makefile-MPICH-ifort-okeanos \ No newline at end of file +Makefile-MPICH-ifort \ No newline at end of file diff --git a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort index 79b8d0f..907ce62 100644 --- a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort +++ b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort @@ -1,10 +1,10 @@ -INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh +INSTALL_DIR = /users/software/mpich2-1.4.1p1_intel BIN=../../../../bin/cluster -FC = ifort +FC= ${INSTALL_DIR}/bin/mpif90 OPT = -O3 -ip -w -mcmodel=medium -OPT = -CB -g -mcmodel=medium +#OPT = -CB -g -mcmodel=medium FFLAGS = ${OPT} -c -I. -Iinclude_unres -I$(INSTALL_DIR)/include -LIBS = -L$(INSTALL_DIR)/lib -lmpich -lpmpich xdrf/libxdrf.a +LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a .c.o: cc -c -DLINUX -DPGI $*.c @@ -20,49 +20,97 @@ object = main_clust.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o geomout.o readpdb.o read_coords.o parmread.o probabl.o fitsq.o hc.o \ track.o wrtclust.o srtclust.o noyes.o contact.o printmat.o \ int_from_cart1.o energy_p_new.o icant.o proc_proc.o work_partition.o \ - setup_var.o read_ref_str.o gnmr1.o permut.o ssMD.o + setup_var.o read_ref_str.o gnmr1.o permut.o seq2chains.o \ + chain_symmetry.o iperm.o rmscalc.o rmsnat.o TMscore.o ssMD.o refsys.o \ + read_constr_homology.o boxshift.o all: no_option - @echo "Specify force field: GAB, 4P or E0LL2Y" + @echo "Specify force field: GAB, 4P, E0LL2Y or NEWCORR" no_option: GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ - -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -GAB: BIN = ../../../../bin/cluster/unres_clustMD-mult_ifort_MPICH_GAB.exe + -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC \ + -DFOURBODY +GAB: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_GAB-HCD.exe GAB: ${object} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c ./compinfo | true ${FC} ${FFLAGS} cinfo.f $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} 4P: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ - -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -4P: BIN = ../../../../bin/cluster/unres_clustMD-mult_ifort_MPICH_4P.exe + -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC \ + -DFOURBODY +4P: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_4P-HCD.exe 4P: ${object} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c ./compinfo | true ${FC} ${FFLAGS} cinfo.f $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} E0LL2Y: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ - -DCLUST -DSPLITELE -DLANG0 -E0LL2Y: BIN = ../../../../bin/cluster/unres_clustMD-mult_ifort_MPICH_E0LL2Y.exe + -DCLUST -DSPLITELE -DFOURBODY +E0LL2Y: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_E0LL2Y-HCD.exe E0LL2Y: ${object} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c ./compinfo | true ${FC} ${FFLAGS} cinfo.f $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} +E0LL2Y_DFA: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCLUST -DSPLITELE -DFOURBODY -DDFA +E0LL2Y_DFA: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_E0LL2Y-HCD-DFA.exe +E0LL2Y_DFA: ${object} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} dfa.o cinfo.o ${LIBS} -o ${BIN} + NEWCORR: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ - -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -NEWCORR: BIN = ../../../../bin/cluster/unres_clustMD-mult_ifort_MPICH_NEWCORR.exe + -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_SC-HCD.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe NEWCORR: ${object} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} + +NEWCORR5D: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -DFIVEDIAG +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR5D: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_SC-HCD5.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe +NEWCORR5D: ${object} xdrf/libxdrf.a + gcc -o compinfo compinfo.c ./compinfo | true ${FC} ${FFLAGS} cinfo.f $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} +NEWCORR_DFA: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -DDFA +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR_DFA: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_SC-HCD-DFA.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe +NEWCORR_DFA: ${object} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} dfa.o cinfo.o ${LIBS} -o ${BIN} + +NEWCORR5D_DFA: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -DDFA +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR5D_DFA: BIN = ~/bin/unres_clustMD_ifort_MPICH-okeanos_SC-HCD5-DFA.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe +NEWCORR5D_DFA: ${object} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} dfa.o cinfo.o ${LIBS} -o ${BIN} + xdrf/libxdrf.a: cd xdrf && make diff --git a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-prometheus b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-prometheus index 1492755..e1b8f32 100644 --- a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-prometheus +++ b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-prometheus @@ -1,14 +1,14 @@ -FC = mpif90 -fc=ifort +################################################################### +#INSTALL_DIR = /net/software/local/intel/compilers_and_libraries_2016.3.210/linux/mpi/intel64 + -OPT = -O3 -ip -mcmodel=medium -shared-intel -#OPT = -O3 -#OPT = -g -CA -CB -mcmodel=medium -shared-intel +FC = mpif90 -fc=ifort -FFLAGS = -c ${OPT} -Iinclude_unres -FFLAGS1 = -c -g -CA -CB -mcmodel=medium -shared-intel -#FFLAGS = ${FFLAGS1} -LIBS = -lmpi xdrf/libxdrf.a +OPT = -O3 -ip -mcmodel=medium +#OPT = -CB -g -mcmodel=medium -shared-intel +FFLAGS = ${OPT} -c -I. -Iinclude_unres -I$(INSTALL_DIR)/include +LIBS = -L$(INSTALL_DIR)/lib -lmpi xdrf/libxdrf.a .c.o: cc -c -DLINUX -DPGI $*.c @@ -24,16 +24,19 @@ object = main_clust.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o geomout.o readpdb.o read_coords.o parmread.o probabl.o fitsq.o hc.o \ track.o wrtclust.o srtclust.o noyes.o contact.o printmat.o \ int_from_cart1.o energy_p_new.o icant.o proc_proc.o work_partition.o \ - setup_var.o read_ref_str.o gnmr1.o permut.o rmsnat.o TMscore.o ssMD.o oligomer.o + setup_var.o read_ref_str.o gnmr1.o permut.o seq2chains.o \ + chain_symmetry.o iperm.o rmscalc.o rmsnat.o TMscore.o ssMD.o refsys.o \ + read_constr_homology.o all: no_option - @echo "Specify force field: GAB, 4P or E0LL2Y" + @echo "Specify force field: GAB, 4P, E0LL2Y or NEWCORR" no_option: GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ - -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -GAB: BIN = ~/unres/bin/unres_clustMD-mult_ifort_MPICH_GAB-SAXS-MRAMB-Bfac.exe + -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC \ + -DFOURBODY +GAB: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_GAB-HCD.exe GAB: ${object} xdrf/libxdrf.a gcc -o compinfo compinfo.c ./compinfo | true @@ -41,8 +44,9 @@ GAB: ${object} xdrf/libxdrf.a $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} 4P: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ - -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -4P: BIN = ~/unres/bin/unres_clustMD-mult_ifort_MPICH_4P-SAXS-MRSAMB-Bfac.exe + -DCLUST -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC \ + -DFOURBODY +4P: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_4P-HCD.exe 4P: ${object} xdrf/libxdrf.a gcc -o compinfo compinfo.c ./compinfo | true @@ -50,23 +54,67 @@ GAB: ${object} xdrf/libxdrf.a $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} E0LL2Y: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ - -DCLUST -DSPLITELE -DLANG0 -E0LL2Y: BIN = ~/unres/bin/unres_clustMD-mult_ifort_MPICH_E0LL2Y-SAXS-MRAMB-Bfac.exe + -DCLUST -DSPLITELE -DFOURBODY +E0LL2Y: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_E0LL2Y-HCD.exe E0LL2Y: ${object} xdrf/libxdrf.a gcc -o compinfo compinfo.c ./compinfo | true ${FC} ${FFLAGS} cinfo.f $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} +E0LL2Y_DFA: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCLUST -DSPLITELE -DFOURBODY -DDFA +E0LL2Y_DFA: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_E0LL2Y-HCD-DFA.exe +E0LL2Y_DFA: ${object} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} dfa.o cinfo.o ${LIBS} -o ${BIN} + NEWCORR: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -NEWCORR: BIN = ~/unres/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-SAXS-MRAMB-Bfac.exe +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_SC-HCD.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe NEWCORR: ${object} xdrf/libxdrf.a gcc -o compinfo compinfo.c ./compinfo | true ${FC} ${FFLAGS} cinfo.f $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} +NEWCORR5D: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -DFIVEDIAG +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR5D: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_SC-HCD5.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe +NEWCORR5D: ${object} xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN} + +NEWCORR_DFA: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -DDFA +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR_DFA: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_SC-HCD-DFA.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe +NEWCORR_DFA: ${object} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} dfa.o cinfo.o ${LIBS} -o ${BIN} + +NEWCORR5D_DFA: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ + -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -DDFA +#-DCLUST -DSPLITELE -DLANG0 -DNEWCORR +NEWCORR5D_DFA: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-prometheus_SC-HCD5-DFA.exe +#NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe +NEWCORR5D_DFA: ${object} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo | true + ${FC} ${FFLAGS} cinfo.f + $(FC) ${OPT} ${object} dfa.o cinfo.o ${LIBS} -o ${BIN} + xdrf/libxdrf.a: cd xdrf && make diff --git a/source/cluster/wham/src-HCD/Makefile-tryton b/source/cluster/wham/src-HCD/Makefile-tryton index e887bc9..ec6ae53 100644 --- a/source/cluster/wham/src-HCD/Makefile-tryton +++ b/source/cluster/wham/src-HCD/Makefile-tryton @@ -85,7 +85,7 @@ NEWCORR: ${object} xdrf/libxdrf.a NEWCORR5D: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI -DPROCOR \ -DCORRCD -DCLUST -DSPLITELE -DLANG0 -DNEWCORR -DFIVEDIAG #-DCLUST -DSPLITELE -DLANG0 -DNEWCORR -NEWCORR5D: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-tryton_SC-HCD5.exe +NEWCORR5D: BIN = ~/unres/bin/unres_clustMD_ifort_MPICH-tryton_SC-HCD5-L.exe #NEWCORR: BIN = ~/bin/unres_clustMD-mult_ifort_MPICH_NEWCORR-fouriertor-test.exe NEWCORR5D: ${object} xdrf/libxdrf.a gcc -o compinfo compinfo.c diff --git a/source/cluster/wham/src-HCD/energy_p_new.F b/source/cluster/wham/src-HCD/energy_p_new.F index db2e043..4fa79c5 100644 --- a/source/cluster/wham/src-HCD/energy_p_new.F +++ b/source/cluster/wham/src-HCD/energy_p_new.F @@ -1104,6 +1104,7 @@ C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include "DIMENSIONS.COMPAR" + include 'COMMON.CONTROL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.LOCAL' @@ -1264,6 +1265,8 @@ C#define DEBUG #endif C#undef DEBUG c endif + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij if (calc_grad) then C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1272,6 +1275,12 @@ C Calculate gradient components. fac=rij*fac fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -1433,6 +1442,12 @@ C Calculate gradient components. fac=rij*fac-2*expon*rrij*e_augm fac=fac+(evdwij+e_augm)*sssgrad/sss*rij C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -2092,6 +2107,8 @@ C common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -2200,6 +2217,7 @@ c end if ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -2230,6 +2248,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2270,6 +2289,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2340,6 +2360,9 @@ C------------------------------------------------------------------------------- common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij, + & faclipij2 + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -2375,6 +2398,9 @@ C zj=c(3,j)+0.5D0*dzj-zmedi yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0 + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -2413,25 +2439,25 @@ C fac_shield(j)=0.6 el1=el1*fac_shield(i)**2*fac_shield(j)**2 el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*faclipij2 else fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*faclipij2 endif - evdw1=evdw1+evdwij*sss + evdw1=evdw1+evdwij*sss*faclipij2 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') - &'evdw1',i,j,evdwij - &,iteli,itelj,aaa,evdw1,sss - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - &fac_shield(i),fac_shield(j) + write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') + & 'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss + write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij, + & faclipij2 endif C @@ -2449,9 +2475,10 @@ C * Radial derivatives. First process both termini of the fragment (i,j) * if (calc_grad) then - ggg(1)=facel*xj - ggg(2)=facel*yj - ggg(3)=facel*zj + aux=(facel*sss+rmij*sssgrad*eesij)*faclipij2 + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j @@ -2549,7 +2576,7 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo if (sss.gt.0.0) then - facvdw=facvdw+sssgrad*rmij*evdwij + facvdw=facvdw+sssgrad*rmij*evdwij*faclipij2 ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj @@ -2579,7 +2606,7 @@ cgrad enddo endif ! calc_grad #else C MARYSIA - facvdw=(ev1+evdwij) + facvdw=(ev1+evdwij)*faclipij2 facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel)*sss @@ -2642,7 +2669,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* - & fac_shield(i)**2*fac_shield(j)**2 + & fac_shield(i)**2*fac_shield(j)**2*sss*faclipij2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -2663,11 +2690,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -2899,14 +2926,14 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij c if (eel_loc_ij.ne.0) c & write (iout,'(a4,2i4,8f9.5)')'chuj', c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) - eel_loc=eel_loc+eel_loc_ij*sss + eel_loc=eel_loc+eel_loc_ij C Now derivative over eel_loc if (calc_grad) then if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. @@ -2963,7 +2990,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -2979,7 +3006,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -2992,7 +3019,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -3004,7 +3031,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -3013,12 +3040,12 @@ C Partial derivatives in virtual-bond dihedral angles gamma & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) aux=eel_loc_ij/sss*sssgrad*rmij ggg(1)=aux*xj @@ -3027,7 +3054,7 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) @@ -3040,22 +3067,27 @@ cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) cgrad enddo cgrad enddo C Remaining derivatives of eello + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij + + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij enddo endif ! calc_grad @@ -3321,6 +3353,8 @@ C Third- and fourth-order contributions from turns common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij j=i+2 c write (iout,*) "eturn3",i,j,j1,j2 a_temp(1,1)=a22 @@ -3358,7 +3392,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) if (energy_dec) write (iout,'(6heturn3,2i5,0pf7.3)') i,i+2, @@ -3368,10 +3402,10 @@ C#ifdef NEWCORR C Derivatives in theta gloc(nphi+i,icg)=gloc(nphi+i,icg) & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C#endif C Derivatives in shield mode @@ -3426,14 +3460,14 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Cartesian derivatives do l=1,3 c ghalf1=0.5d0*agg(l,1) @@ -3447,7 +3481,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) @@ -3456,7 +3490,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -3464,7 +3498,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -3472,7 +3506,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif ! calc_grad @@ -3605,7 +3639,7 @@ C fac_shield(i)=0.6 C fac_shield(j)=0.4 endif eello_turn4=eello_turn4-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) @@ -3682,7 +3716,7 @@ C Derivatives in gamma(i) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -3691,7 +3725,7 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -3703,7 +3737,7 @@ C Derivatives in gamma(i+2) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij if (calc_grad) then C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) @@ -3724,7 +3758,7 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif C Remaining derivatives of this turn contribution @@ -3743,7 +3777,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -3758,7 +3792,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -3773,7 +3807,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -3789,7 +3823,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif ! calc_grad diff --git a/source/cluster/wham/src-HCD/include_unres/COMMON.INTERACT b/source/cluster/wham/src-HCD/include_unres/COMMON.INTERACT index 1c0b8db..a02f7e4 100644 --- a/source/cluster/wham/src-HCD/include_unres/COMMON.INTERACT +++ b/source/cluster/wham/src-HCD/include_unres/COMMON.INTERACT @@ -31,6 +31,7 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters. & distchainmax,nbondterm(ntyp) &,vbldpDUM C 01/29/15 Lipidic parameters - double precision pepliptran,liptranene - common /lipid/ pepliptran,liptranene(ntyp) - + double precision pepliptran,liptranene, lipscale, + &tubetranene, tubetranenepep + common /lipid/ pepliptran,liptranene(ntyp),lipscale + common /tubepar/ tubetranene(ntyp), tubetranenepep diff --git a/source/cluster/wham/src-HCD/initialize.f b/source/cluster/wham/src-HCD/initialize.f deleted file mode 100644 index 12ea156..0000000 --- a/source/cluster/wham/src-HCD/initialize.f +++ /dev/null @@ -1,99 +0,0 @@ - subroutine initialize -C -C Define constants and zero out tables. -C - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.IOUNITS' - include 'COMMON.CHAIN' - include 'COMMON.INTERACT' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.TORSION' - include 'COMMON.FFIELD' - include 'COMMON.SBRIDGE' - include 'COMMON.MINIM' - include 'COMMON.DERIV' -C -C The following is just to define auxiliary variables used in angle conversion -C - pi=4.0D0*datan(1.0D0) - dwapi=2.0D0*pi - dwapi3=pi/3.0D0 - pipol=0.5D0*pi - deg2rad=pi/180.0D0 - rad2deg=1.0D0/deg2rad - angmin=10.0D0*deg2rad -C Assign virtual-bond length - vbl=3.8D0 - vblinv=1.0D0/vbl - vblinv2=vblinv*vblinv -C -C Define I/O units. -C - inp= 1 - iout= 2 - ipdbin= 3 - ipdb= 7 - igeom= 8 - intin= 9 - istat= 17 - imol2= 18 - jplot= 19 - jstatin=10 - jstatout=11 -C -C Zero out tables. -C - do i=1,maxres2 - do j=1,3 - c(j,i)=0.0D0 - dc(j,i)=0.0D0 - enddo - enddo - do i=1,maxres - do j=1,3 - xloc(j,i)=0.0D0 - enddo - enddo -C Initialize the bridge arrays - ns=0 - nss=0 - nhpb=0 - do i=1,maxss - iss(i)=0 - enddo - do i=1,maxdim - dhpb(i)=0.0D0 - enddo - do i=1,maxres - ihpb(i)=0 - jhpb(i)=0 - enddo -C -C Initialize timing. -C - call set_timers - return - end -c------------------------------------------------------------------------- - block data chuj - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.NAMES' - include 'COMMON.FFIELD' - data restyp / - &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL', - & 'DSG','DGN','DSN','DTH', - &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER', - &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR', - &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ', - &'AIB','ABU','D'/ - data onelet / - &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g', - &'a','y','w','v','l','i','f','m','c','x', - &'C','M','F','I','L','V','W','Y','A','G','T', - &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/ - data potname /'LJ','LJK','BP','GB','GBV'/ - data potname /'LJ','LJK','BP','GB','GBV'/ - end diff --git a/source/cluster/wham/src-HCD/readrtns.F b/source/cluster/wham/src-HCD/readrtns.F index 057f1ac..e9e576f 100644 --- a/source/cluster/wham/src-HCD/readrtns.F +++ b/source/cluster/wham/src-HCD/readrtns.F @@ -234,6 +234,7 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'WDFAN',wdfa_nei,0.0d0) call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) call reada(weightcard,'WLT',wliptran,0.0D0) + call reada(weightcard,'LIPSCALE',lipscale,1.0D0) call reada(weightcard,"ATRISS",atriss,0.301D0) call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) diff --git a/source/cluster/wham/src-HCD/sizesclu.dat b/source/cluster/wham/src-HCD/sizesclu.dat index cb8572c..7d0d666 100644 --- a/source/cluster/wham/src-HCD/sizesclu.dat +++ b/source/cluster/wham/src-HCD/sizesclu.dat @@ -5,7 +5,7 @@ * Max. number of conformations in the data set. * integer maxconf,maxstr_proc - PARAMETER (MAXCONF=12000) + PARAMETER (MAXCONF=8000) parameter (maxstr_proc=maxconf/2) * * Max. number of "distances" between conformations. diff --git a/source/unres/src-HCD-5D/COMMON.DFA b/source/unres/src-HCD-5D/COMMON.DFA index 67a1a0d..6759f8b 100644 --- a/source/unres/src-HCD-5D/COMMON.DFA +++ b/source/unres/src-HCD-5D/COMMON.DFA @@ -68,14 +68,21 @@ C NMAP - mapping between dfanum and ndis, nphi, nthe, nnei & NCA, ICAIDX(MAXRES) COMMON /IDFA2/ ishiftca,ilastca c parallel - integer idfadis_start,idfadis_end, - & idfaphi_start,idfaphi_end, - & idfathe_start,idfathe_end, - & idfanei_start,idfanei_end + integer idfadis_start,idfadis_end,idfaphi_start,idfaphi_end, + & idfathe_start,idfathe_end,idfanei_start,idfanei_end, + & idfadis_start_all(0:max_fg_procs-1), + & idfadis_end_all(0:max_fg_procs-1), + & idfaphi_start_all(0:max_fg_procs-1), + & idfaphi_end_all(0:max_fg_procs-1), + & idfathe_start_all(0:max_fg_procs-1), + & idfathe_end_all(0:max_fg_procs-1), + & idfanei_start_all(0:max_fg_procs-1), + & idfanei_end_all(0:max_fg_procs-1) COMMON /dfa_mpi/ idfadis_start,idfadis_end, - & idfaphi_start,idfaphi_end, - & idfathe_start,idfathe_end, - & idfanei_start,idfanei_end + & idfaphi_start,idfaphi_end,idfathe_start,idfathe_end, + & idfanei_start,idfanei_end,idfadis_start_all,idfadis_end_all, + & idfaphi_start_all,idfaphi_end_all,idfathe_start_all, + & idfathe_end_all,idfanei_start_all,idfanei_end_all CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC diff --git a/source/unres/src-HCD-5D/COMMON.INTERACT b/source/unres/src-HCD-5D/COMMON.INTERACT index 8c4876d..14416ad 100644 --- a/source/unres/src-HCD-5D/COMMON.INTERACT +++ b/source/unres/src-HCD-5D/COMMON.INTERACT @@ -58,8 +58,8 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters. & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp), & distchainmax,nbondterm(ntyp) C 01/29/15 Lipidic parameters - double precision pepliptran,liptranene, + double precision pepliptran,liptranene, lipscale, &tubetranene, tubetranenepep - common /lipid/ pepliptran,liptranene(ntyp) + common /lipid/ pepliptran,liptranene(ntyp),lipscale common /tubepar/ tubetranene(ntyp), tubetranenepep diff --git a/source/unres/src-HCD-5D/DIMENSIONS b/source/unres/src-HCD-5D/DIMENSIONS index 137b45d..9803b23 100644 --- a/source/unres/src-HCD-5D/DIMENSIONS +++ b/source/unres/src-HCD-5D/DIMENSIONS @@ -16,10 +16,10 @@ C Max. number of coarse-grain processors parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres - parameter (maxres=10000) + parameter (maxres=2000) C Max. number of AA residues per chain integer maxres_chain - parameter (maxres_chain=1200) + parameter (maxres_chain=800) C Max. number of cysteines and other bridging residues integer max_cyst parameter (max_cyst=100) diff --git a/source/unres/src-HCD-5D/MD_A-MTS.F b/source/unres/src-HCD-5D/MD_A-MTS.F index fcef69e..e504cbd 100644 --- a/source/unres/src-HCD-5D/MD_A-MTS.F +++ b/source/unres/src-HCD-5D/MD_A-MTS.F @@ -1955,6 +1955,7 @@ c 8/22/17 AL Loop to produce a low-energy random conformation enddo call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) + dc(:,0)=c(:,1) do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) diff --git a/source/unres/src-HCD-5D/checkder_p.F b/source/unres/src-HCD-5D/checkder_p.F index 48eedda..e1448db 100644 --- a/source/unres/src-HCD-5D/checkder_p.F +++ b/source/unres/src-HCD-5D/checkder_p.F @@ -132,6 +132,7 @@ 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 + call cartprint nf=0 icall=0 call geom_to_var(nvar,x) @@ -212,6 +213,14 @@ c call flush(iout) enddo enddo endif +c write (iout,*) "Vector dc" +c do i=0,nres +c write (iout,'(i5,2(3f10.5,5x))') +c & i,(dc(j,i),j=1,3),(dc(j,i+nres),j=1,3) +c enddo +c write (iout,*) "Coordinates after chainbuild_cart" +c call chainbuild_cart +c call cartprint write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors' do i=0,nres print *,i diff --git a/source/unres/src-HCD-5D/dfa.F b/source/unres/src-HCD-5D/dfa.F index f69b81a..62e8892 100644 --- a/source/unres/src-HCD-5D/dfa.F +++ b/source/unres/src-HCD-5D/dfa.F @@ -71,12 +71,16 @@ C read fragment informations C implicit real*8 (a-h,o-z) include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" + include 'COMMON.SETUP' + integer ierror +#endif include 'COMMON.IOUNITS' include 'COMMON.CHAIN' include 'COMMON.DFA' include 'COMMON.FFIELD' include 'COMMON.CONTROL' - include 'COMMON.SETUP' C NOTE THAT FILENAMES are FIXED, CURRENTLY!! @@ -90,6 +94,7 @@ C THIS SHOULD BE MODIFIED!! integer ica1, ica2,ica3,ica4,ica5 integer ishell, inca, itmp,iitmp double precision wtmp + logical lprn /.false./ C C READ DISTANCE C @@ -254,12 +259,55 @@ C BETA is not parallel ! call int_bounds(idfaphi,idfaphi_start,idfaphi_end) call int_bounds(idfathe,idfathe_start,idfathe_end) call int_bounds(idfanei,idfanei_start,idfanei_end) - if (me.eq.king .or. .not. out1file) + if (lprn) write (*,*) "Processor",MyRank," DFA MPI ", + & "idfadis ",idfadis,idfadis_start,idfadis_end, + & "idfaphi ",idfaphi,idfaphi_start,idfaphi_end, + & "idfathe ",idfathe,idfathe_start,idfathe_end, + & "idfanei ",idfanei,idfanei_start,idfanei_end + if (lprn) & write (iout,*) "DFA MPI ", & "idfadis ",idfadis,idfadis_start,idfadis_end, & "idfaphi ",idfaphi,idfaphi_start,idfaphi_end, & "idfathe ",idfathe,idfathe_start,idfathe_end, & "idfanei ",idfanei,idfanei_start,idfanei_end + do i=0,max_fg_procs-1 + idfadis_start_all(j)=0 + idfadis_end_all(j)=0 + idfaphi_start_all(j)=0 + idfaphi_end_all(j)=0 + idfathe_start_all(j)=0 + idfathe_end_all(j)=0 + idfanei_start_all(j)=0 + idfanei_end_all(j)=0 + enddo + call MPI_Allgather(idfadis_start,1,MPI_INTEGER, + & idfadis_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgather(idfadis_end,1,MPI_INTEGER, + & idfadis_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgather(idfaphi_start,1,MPI_INTEGER, + & idfaphi_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgather(idfaphi_end,1,MPI_INTEGER, + & idfaphi_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgather(idfathe_start,1,MPI_INTEGER, + & idfathe_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgather(idfathe_end,1,MPI_INTEGER, + & idfathe_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgather(idfanei_start,1,MPI_INTEGER, + & idfanei_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgather(idfanei_end,1,MPI_INTEGER, + & idfanei_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR) + if (me.eq.0 .or. out1file) then + write (iout,*) "Partitioning of DFA work" + write (iout,'(5a10)') 'Rank','DFA_dis','DFA_phi','DFA_the', + & 'DFA_nei' + do i=0,nfgtasks-1 + write (iout,'(i10,8i5)') i,idfadis_start_all(i), + & idfadis_end_all(i),idfaphi_start_all(i), + & idfaphi_end_all(i),idfathe_start_all(i), + & idfathe_end_all(i),idfanei_start_all(i), + & idfanei_end_all(i) + enddo + endif #else idfadis_start=1 idfadis_end=idfadis diff --git a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F index 28ba1d1..c4e54bc 100644 --- a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F @@ -84,10 +84,15 @@ c include 'COMMON.CONTACTS' integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sss1,sssgrad1,sqrij - double precision sscale,sscagrad + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision sscale,sscagrad,sscagradlip,sscalelip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -99,6 +104,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -112,6 +118,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -127,6 +138,7 @@ c do j=istart(i,iint),iend(i,iint) if (sss.lt.1.0d0) then rrij=1.0D0/rij fac=rrij**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e1+e2 @@ -140,11 +152,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -192,10 +209,15 @@ c include 'COMMON.CONTACTS' integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -207,6 +229,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C Change 12/1/95 num_conti=0 C @@ -222,6 +245,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -235,6 +263,7 @@ C Change 12/1/95 to calculate four-body interactions rrij=1.0D0/rij eps0ij=eps(itypi,itypj) fac=rrij**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e1+e2 @@ -246,11 +275,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -296,10 +330,15 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -311,6 +350,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -322,6 +362,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -339,6 +384,7 @@ c do j=istart(i,iint),iend(i,iint) & sscagrad(rij/sigma(itypi,itypj),r_cut_respa) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e_augm+e1+e2 @@ -360,11 +406,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -401,10 +452,15 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift + double precision gg_lipi(3),gg_lipj(3) c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -416,6 +472,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -427,6 +484,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -439,6 +501,7 @@ c do j=istart(i,iint),iend(i,iint) if (sss.gt.0.0d0) then r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=e_augm+e1+e2 @@ -459,11 +522,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo endif c enddo ! j @@ -501,13 +569,16 @@ C integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision sss1,sssgrad1 - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c if (icall.eq.0) then c lprn=.true. c else @@ -525,6 +596,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -553,6 +625,11 @@ c dscj_inv=dsc_inv(itypj) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -573,6 +650,7 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives fac=(rrij*sigsq)**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -600,6 +678,12 @@ C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad_scale((1.0d0-sss)*sss1) @@ -633,11 +717,15 @@ C double precision evdw integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi - double precision sscale,sscagrad + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) then c lprn=.true. @@ -656,6 +744,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -684,6 +773,11 @@ c dscj_inv=dsc_inv(itypj) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -701,6 +795,7 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives fac=(rrij*sigsq)**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -727,6 +822,17 @@ C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi + do k=1,3 + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) + enddo C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad_scale(sss) @@ -770,7 +876,8 @@ C evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 @@ -822,9 +929,9 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) call to_box(xj,yj,zj) call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -860,6 +967,7 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -880,8 +988,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,4f10.5)') - & 'evdw',i,j,rij,sss,sss1,evdwij + if (energy_dec) write (iout,'(a,2i5,5f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss1,sss,sslipi,sslipj,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -894,8 +1002,13 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac - gg_lipi(3)=ssgradlipi*evdwij - gg_lipj(3)=ssgradlipj*evdwij + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi + C Calculate angular part of the gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif @@ -936,7 +1049,8 @@ C evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0D0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 @@ -988,9 +1102,14 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) call to_box(xj,yj,zj) call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 +c write (iout,*) "aa bb",aa_lip(itypi,itypj), +c & bb_lip(itypi,itypj),aa_aq(itypi,itypj), +c & bb_aq(itypi,itypj),aa,bb +c write (iout,*) (sslipi+sslipj)/2.0d0, +c & (2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1023,6 +1142,7 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1043,8 +1163,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1056,8 +1176,13 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac - gg_lipi(3)=ssgradlipi*evdwij - gg_lipj(3)=ssgradlipj*evdwij + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi +c write (iout,*) "gglip",i,j,gg_lipi,gg_lipj C Calculate angular part of the gradient. call sc_grad_scale(sss) endif @@ -1098,6 +1223,8 @@ C double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision sss1,sssgrad1 evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.true. @@ -1181,6 +1308,7 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1213,6 +1341,12 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate angular part of the gradient. call sc_grad_scale((1.0d0-sss)*sss1) endif @@ -1251,6 +1385,8 @@ C double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision boxshift evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.true. @@ -1330,6 +1466,7 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1360,6 +1497,12 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate angular part of the gradient. call sc_grad_scale(sss) endif @@ -1414,6 +1557,7 @@ c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv C C Calculate the components of the gradient in DC and X C +c write (iout,*) "scgrad gglip",i,j,gg_lipi,gg_lipj do l=1,3 gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) @@ -1462,6 +1606,8 @@ C common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -1556,6 +1702,7 @@ C & .or. itype(i+4).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 call eelecij_scale(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -1580,6 +1727,7 @@ C & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -1611,6 +1759,7 @@ C & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend #ifdef FOURBODY num_conti=num_cont_hb(i) @@ -1693,7 +1842,9 @@ C------------------------------------------------------------------------------- double precision sscale,sscagrad double precision scalar double precision boxshift - + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij, + & faclipij2 + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -1726,6 +1877,9 @@ C print *,"WCHODZE2" yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0 + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -1763,17 +1917,18 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions eesij=el1+el2 C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) - ees=ees+eesij*sss1 - evdw1=evdw1+evdwij*(1.0d0-sss)*sss1 + ees=ees+eesij*sss1*faclipij2 + evdw1=evdw1+evdwij*(1.0d0-sss)*sss1*faclipij2 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2f7.3)') - & 'evdw1',i,j,evdwij,sss,sss1 - write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij + write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,5f10.5)') + & 'evdw1',i,j,evdwij,iteli,itelj,aaa,sss,sss1,sssgrad,sssgrad1,rij + write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij,faclipij2 endif C @@ -1791,7 +1946,8 @@ c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) * * Radial derivatives. First process both termini of the fragment (i,j) * - aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij +c old aux=(facel+sssgrad1*(1.0d0-sss)*eesij*rmij)*faclipij2 + aux=(facel+sssgrad1*eesij*rmij)*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(1)=aux*xj ggg(2)=aux*yj @@ -1864,6 +2020,10 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2*sss1 + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2*sss1 c gelc_long(3,i)=gelc_long(3,i)+ c ssgradlipi*eesij/2.0d0*lipscale**2*sss1 * @@ -1874,9 +2034,9 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - facvdw=facvdw+ - & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij -c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + facvdw=(facvdw+ + &(-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij) + & *faclipij2 ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj @@ -1890,6 +2050,11 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss1*(1.0d0-sss)*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss1*(1.0d0-sss)*ssgradlipi*evdwij/2.0d0*lipscale**2 * * Loop over residues i+1 thru j-1. * @@ -1914,8 +2079,8 @@ c facel=el1+eesij * * Radial derivatives. First process both termini of the fragment (i,j) * - aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj)) - & *eesij*rmij + aux=(fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj)) + & *eesij*rmij)*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) ggg(1)=aux*xj ggg(2)=aux*yj @@ -1954,6 +2119,10 @@ C ggg(3)=facvdw*zj gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss1*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss1*ssgradlipi*evdwij/2.0d0*lipscale**2 #endif * * Angular part @@ -1971,7 +2140,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1 - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) enddo @@ -1993,13 +2162,13 @@ cgrad enddo gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1 - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & +ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss1 - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) @@ -2213,7 +2382,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij eel_loc=eel_loc+eel_loc_ij if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. @@ -2265,7 +2434,7 @@ C & *2.0 & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -2281,7 +2450,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -2294,7 +2463,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -2306,14 +2475,14 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij #endif cC Paral derivatives in virtual-bond dihedral angles gamma if (i.gt.1) & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) @@ -2321,7 +2490,7 @@ c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) @@ -2333,7 +2502,7 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *fac_shield(i)*fac_shield(j) c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) @@ -2343,6 +2512,10 @@ cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -2360,22 +2533,22 @@ c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss1 + & *fac_shield(i)*fac_shield(j)*sss1*faclipij c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo @@ -2629,8 +2802,10 @@ c write (iout,*) "evdwpp_short" double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, & dist_temp, dist_init,sss_grad double precision sscale,sscagrad + double precision sslipi,ssgradlipi,sslipj,ssgradlipj double precision boxshift integer ikont + double precision faclipij2 evdw1=0.0D0 C print *,"WCHODZE" c write (iout,*) "iatel_s_vdw",iatel_s_vdw, @@ -2651,6 +2826,7 @@ c do i=iatel_s_vdw,iatel_e_vdw ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i), c & ' ielend',ielend_vdw(i) @@ -2673,6 +2849,8 @@ c do j=ielstart_vdw(i),ielend_vdw(i) yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -2695,16 +2873,17 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions if (energy_dec) then write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss endif - evdw1=evdw1+evdwij*sss + evdw1=evdw1+evdwij*sss*faclipij2 if (energy_dec) write (iout,'(a10,2i5,0pf7.3)') & 'evdw1_sum',i,j,evdw1 C C Calculate contributions to the Cartesian gradient. C - facvdw=-6*rrmij*(ev1+evdwij)*sss - ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) - ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) - ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) + facvdw=(-6*rrmij*(ev1+evdwij)*sss+sssgrad*rmij*evdwij/ + & rpp(iteli,itelj))*faclipij2 + ggg(1)=facvdw*xj + ggg(2)=facvdw*yj + ggg(3)=facvdw*zj C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj @@ -2712,6 +2891,11 @@ C ggg(3)=facvdw*zj gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 endif c enddo ! j enddo ! i diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F index 2c701ca..6d5b25f 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -103,10 +103,10 @@ C FG slaves receive the WEIGHTS array wliptran=weights(22) wtube=weights(25) wsaxs=weights(26) - wdfa_dist=weights_(28) - wdfa_tor=weights_(29) - wdfa_nei=weights_(30) - wdfa_beta=weights_(31) + wdfa_dist=weights(28) + wdfa_tor=weights(29) + wdfa_nei=weights(30) + wdfa_beta=weights(31) endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 @@ -177,8 +177,10 @@ C 107 continue #ifdef DFA C BARTEK for dfa test! +c print *,"Processors",MyRank," wdfa",wdfa_dist if (wdfa_dist.gt.0) then call edfad(edfadis) +c print *,"Processors",MyRank," edfadis",edfadis else edfadis=0 endif @@ -831,7 +833,8 @@ c call flush(iout) #ifdef TIMING c time_allreduce=time_allreduce+MPI_Wtime()-time00 #endif - do i=nnt,nres +c do i=nnt,nres + do i=0,nres do k=1,3 gradbufc(k,i)=0.0d0 enddo @@ -856,7 +859,8 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,-1,-1 +c do i=nres-2,-1,-1 + do i=nres-2,0,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -872,12 +876,13 @@ c enddo #endif #ifdef DEBUG write (iout,*) "gradbufc" - do i=1,nres + do i=0,nres write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) enddo call flush(iout) #endif - do i=-1,nres +c do i=-1,nres + do i=0,nres do j=1,3 gradbufc_sum(j,i)=gradbufc(j,i) gradbufc(j,i)=0.0d0 @@ -886,7 +891,8 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,-1,-1 +c do i=nres-2,-1,-1 + do i=nres-2,0,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -914,7 +920,8 @@ c enddo do k=1,3 gradbufc(k,nres)=0.0d0 enddo - do i=-1,nct +c do i=-1,nct + do i=0,nct do j=1,3 #ifdef SPLITELE C print *,gradbufc(1,13) @@ -1019,6 +1026,8 @@ C print *,gradafm(1,13),"AFM" endif #ifdef DEBUG write (iout,*) "gradc gradx gloc after adding" + write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)') + & i,(gradc(j,0,icg),j=1,3),(gradx(j,0,icg),j=1,3) do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)') & i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg) @@ -1048,7 +1057,7 @@ C print *,gradafm(1,13),"AFM" #ifdef MPI if (nfgtasks.gt.1) then do j=1,3 - do i=1,nres + do i=0,nres gradbufc(j,i)=gradc(j,i,icg) gradbufx(j,i)=gradx(j,i,icg) enddo @@ -1075,9 +1084,9 @@ c#undef DEBUG call MPI_Barrier(FG_COMM,IERR) time_barrier_g=time_barrier_g+MPI_Wtime()-time00 time00=MPI_Wtime() - call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres, + call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*(nres+1), & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) - call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres, + call MPI_Reduce(gradbufx(1,0),gradx(1,0,icg),3*(nres+1), & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) @@ -1087,7 +1096,7 @@ c#undef DEBUG time_reduce=time_reduce+MPI_Wtime()-time00 #ifdef DEBUG write (iout,*) "gradc after reduce" - do i=1,nres + do i=0,nres do j=1,3 write (iout,*) i,j,gradc(j,i,icg) enddo @@ -1496,10 +1505,15 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 double precision fcont,fprimcont - double precision sscale,sscagrad + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision sscale,sscagrad,sscagradlip,sscalelip + double precision gg_lipi(3),gg_lipj(3) double precision boxshift c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -1511,6 +1525,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C Change 12/1/95 num_conti=0 C @@ -1526,6 +1541,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1540,6 +1560,7 @@ C Change 12/1/95 to calculate four-body interactions c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 + faclip=fac C have you changed here? e1=fac*fac*aa e2=fac*bb @@ -1559,11 +1580,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo cgrad do k=i,j-1 cgrad do l=1,3 @@ -1675,10 +1701,15 @@ C double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck - double precision sscale,sscagrad double precision boxshift + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision gg_lipi(3),gg_lipj(3) + double precision sscale,sscagrad,sscagradlip,sscalelip c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e do ikont=g_listscsc_start,g_listscsc_end i=newcontlisti(ikont) @@ -1690,6 +1721,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) C C Calculate SC interaction energy. C @@ -1701,6 +1733,11 @@ c do j=istart(i,iint),iend(i,iint) yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1714,6 +1751,7 @@ c do j=istart(i,iint),iend(i,iint) sssgrad1=sscagrad(rij,r_cut_int) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon + faclip=fac C have you changed here? e1=fac*fac*aa e2=fac*bb @@ -1734,11 +1772,16 @@ C gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=(sss1/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) - gvdwx(k,j)=gvdwx(k,j)+gg(k) - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k) enddo cgrad do k=i,j-1 cgrad do l=1,3 @@ -1780,13 +1823,16 @@ C integer itypi,itypj,itypi1,iint,ind,ikont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi, & sss1,sssgrad1 - double precision sscale,sscagrad + double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi, + & faclip + double precision sscale,sscagrad,sscagradlip,sscalelip double precision boxshift c double precision rrsave(maxdim) logical lprn evdw=0.0D0 c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 c if (icall.eq.0) then c lprn=.true. c else @@ -1804,6 +1850,7 @@ c do i=iatsc_s,iatsc_e yi=c(2,nres+i) zi=c(3,nres+i) call to_box(xi,yi,zi) + call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1842,6 +1889,11 @@ c alf12=0.0D0 yj=c(2,nres+j) zj=c(3,nres+j) call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + 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 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1864,6 +1916,7 @@ C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives C have you changed here? fac=(rrij*sigsq)**expon2 + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1891,6 +1944,12 @@ C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss1/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad @@ -1931,7 +1990,8 @@ C evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 @@ -2043,9 +2103,15 @@ c alf12=0.0D0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 -C write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) -C if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)') -C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) +c write (iout,*) "aa bb",aa_lip(itypi,itypj), +c & bb_lip(itypi,itypj),aa_aq(itypi,itypj), +c & bb_aq(itypi,itypj),aa,bb +c write (iout,*) (sslipi+sslipj)/2.0d0, +c & (2.0d0-sslipi-sslipj)/2.0d0 + +c write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) +c if (aa.ne.aa_aq(itypi,itypj)) write(iout,'(2e15.5)') +c &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) C if (ssgradlipj.gt.0.0d0) print *,"??WTF??" C print *,sslipi,sslipj,bordlipbot,zi,zj xj=boxshift(xj-xi,boxxsize) @@ -2116,8 +2182,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a,2i5,2f10.5,e15.5)') - & 'r sss evdw',i,j,1.0d0/rij,sss,evdwij + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -2185,7 +2251,8 @@ C double precision dist,sscale,sscagrad,sscagradlip,sscalelip evdw=0.0D0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon - evdw=0.0D0 + gg_lipi=0.0d0 + gg_lipj=0.0d0 lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 @@ -2281,6 +2348,7 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -2310,7 +2378,7 @@ C Calculate gradient components. C Calculate the radial part of the gradient gg_lipi(3)=eps1*(eps2rt*eps2rt) & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* - & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) gg_lipj(3)=ssgradlipj*gg_lipi(3) gg_lipi(3)=gg_lipi(3)*ssgradlipi @@ -3479,6 +3547,8 @@ C common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -3585,6 +3655,7 @@ c end if ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -3640,6 +3711,7 @@ c & (zmedi.lt.((-0.5d0)*boxzsize))) then c go to 196 c endif call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -3683,6 +3755,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) C xmedi=xmedi+xshift*boxxsize C ymedi=ymedi+yshift*boxysize C zmedi=zmedi+zshift*boxzsize @@ -3802,6 +3875,9 @@ C------------------------------------------------------------------------------- double precision xmedi,ymedi,zmedi double precision sscale,sscagrad,scalar double precision boxshift + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij, + & faclipij2 + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -3816,6 +3892,7 @@ C 13-go grudnia roku pamietnego... c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j c ind=ind+1 +c write (iout,*) "lipscale",lipscale iteli=itel(i) itelj=itel(j) if (j.eq.i+2 .and. itelj.eq.2) iteli=2 @@ -3836,6 +3913,9 @@ C zj=c(3,j)+0.5D0*dzj-zmedi yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0 + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -3875,14 +3955,15 @@ C fac_shield(j)=0.6 el1=el1*fac_shield(i)**2*fac_shield(j)**2 el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*sss*faclipij2 else fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) - ees=ees+eesij*sss + ees=ees+eesij*sss*faclipij2 endif - evdw1=evdw1+evdwij*sss + ees=ees + evdw1=evdw1+evdwij*sss*faclipij2 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, @@ -3891,8 +3972,9 @@ cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,3f10.5)') & 'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss,rij - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - & fac_shield(i),fac_shield(j) + write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij, + & faclipij2 endif C @@ -3909,7 +3991,7 @@ C * * Radial derivatives. First process both termini of the fragment (i,j) * - aux=facel*sss+rmij*sssgrad*eesij + aux=(facel*sss+rmij*sssgrad*eesij)*faclipij2 ggg(1)=aux*xj ggg(2)=aux*yj ggg(3)=aux*zj @@ -3991,15 +4073,14 @@ c 9/28/08 AL Gradient compotents will be summed only at the end C print *,"before", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc_long(k,j)=gelc_long(k,j)+ggg(k) -C & +grad_shield(k,j)*eesij/fac_shield(j) gelc_long(k,i)=gelc_long(k,i)-ggg(k) -C & +grad_shield(k,i)*eesij/fac_shield(i) -C gelc_long(k,i-1)=gelc_long(k,i-1) -C & +grad_shield(k,i)*eesij/fac_shield(i) -C gelc_long(k,j-1)=gelc_long(k,j-1) -C & +grad_shield(k,j)*eesij/fac_shield(j) enddo -C print *,"bafter", gelc_long(1,i), gelc_long(1,j) + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2*sss + + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2*sss + * * Loop over residues i+1 thru j-1. @@ -4009,7 +4090,7 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - facvdw=facvdw+sssgrad*rmij*evdwij + facvdw=(facvdw+sssgrad*rmij*evdwij)*faclipij2 ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj @@ -4023,6 +4104,11 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 * * Loop over residues i+1 thru j-1. * @@ -4033,7 +4119,7 @@ cgrad enddo cgrad enddo #else C MARYSIA - facvdw=(ev1+evdwij) + facvdw=(ev1+evdwij)*faclipij2 facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel)*sss @@ -4076,6 +4162,10 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 #endif * * Angular part @@ -4093,7 +4183,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* - & fac_shield(i)**2*fac_shield(j)**2*sss + & fac_shield(i)**2*fac_shield(j)**2*sss*faclipij2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -4114,11 +4204,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -4353,7 +4443,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') c & 'eelloc',i,j,eel_loc_ij C Now derivative over eel_loc @@ -4411,7 +4501,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -4427,7 +4517,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -4440,7 +4530,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -4452,7 +4542,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -4468,12 +4558,12 @@ C Partial derivatives in virtual-bond dihedral angles gamma & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) aux=eel_loc_ij/sss*sssgrad*rmij ggg(1)=aux*xj @@ -4482,13 +4572,19 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij + + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij + cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -4498,19 +4594,19 @@ C Remaining derivatives of eello do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ - & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ - & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ - & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) + & *fac_shield(i)*fac_shield(j)*sss*faclipij enddo ENDIF @@ -4763,6 +4859,8 @@ C Third- and fourth-order contributions from turns common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij j=i+2 c write (iout,*) "eturn3",i,j,j1,j2 a_temp(1,1)=a22 @@ -4800,7 +4898,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) if (energy_dec) write (iout,'(6heturn3,2i5,0pf7.3)') i,i+2, @@ -4809,10 +4907,10 @@ C#ifdef NEWCORR C Derivatives in theta gloc(nphi+i,icg)=gloc(nphi+i,icg) & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C#endif C Derivatives in shield mode @@ -4867,14 +4965,14 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Cartesian derivatives do l=1,3 c ghalf1=0.5d0*agg(l,1) @@ -4888,7 +4986,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) @@ -4897,7 +4995,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -4905,7 +5003,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -4913,8 +5011,17 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo + gshieldc_t3(3,i)=gshieldc_t3(3,i)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j)=gshieldc_t3(3,j)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + return end C------------------------------------------------------------------------------- @@ -5043,7 +5150,7 @@ C fac_shield(i)=0.6 C fac_shield(j)=0.4 endif eello_turn4=eello_turn4-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) @@ -5091,12 +5198,6 @@ C & *2.0 & grad_shield(k,j)*eello_t4/fac_shield(j) enddo endif - - - - - - cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), cd & ' eello_turn4_num',8*eello_turn4_num #ifdef NEWCORR @@ -5126,7 +5227,7 @@ C Derivatives in gamma(i) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -5135,7 +5236,7 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -5147,7 +5248,7 @@ C Derivatives in gamma(i+2) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -5167,7 +5268,7 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif C Remaining derivatives of this turn contribution @@ -5186,7 +5287,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -5201,7 +5302,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -5216,7 +5317,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -5232,8 +5333,16 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) - enddo + & *fac_shield(i)*fac_shield(j)*faclipij + enddo + gshieldc_t4(3,i)=gshieldc_t4(3,i)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j)=gshieldc_t4(3,j)+ + & ssgradlipj*eello_t4/4.0d0*lipscale + gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+ + & ssgradlipj*eello_t4/4.0d0*lipscale return end C----------------------------------------------------------------------------- @@ -11811,6 +11920,7 @@ C--bufliptop--- here true lipid starts C lipid C--buflipbot--- lipid ends buffore starts C--bordlipbot--buffore ends +c call cartprint eliptran=0.0 do i=ilip_start,ilip_end C do i=1,1 @@ -11865,6 +11975,8 @@ CV do i=1,1 if (itype(i).eq.ntyp1) cycle positi=(mod(c(3,i+nres),boxzsize)) if (positi.le.0) positi=positi+boxzsize +c write(iout,*) "i",i," positi",positi,bordlipbot,buflipbot, +c & bordliptop C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop c for each residue check if it is in lipid or lipid water border area C respos=mod(c(3,i+nres),boxzsize) @@ -11875,6 +11987,8 @@ C the energy transfer exist if (positi.lt.buflipbot) then fracinbuf=1.0d0- & ((positi-bordlipbot)/lipbufthick) +c write (iout,*) "i",i,itype(i)," fracinbuf",fracinbuf +c write (iout,*) "i",i," liptranene",liptranene(itype(i)) C lipbufthick is thickenes of lipid buffore sslip=sscalelip(fracinbuf) ssgradlip=-sscagradlip(fracinbuf)/lipbufthick @@ -13250,11 +13364,16 @@ c-------------------------------------------------------------------------- subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi) implicit none include 'DIMENSIONS' + include 'COMMON.IOUNITS' include 'COMMON.CHAIN' double precision xi,yi,zi,sslipi,ssgradlipi double precision fracinbuf double precision sscalelip,sscagradlip - +#ifdef DEBUG + write (iout,*) "bordlipbot",bordlipbot," bordliptop",bordliptop + write (iout,*) "buflipbot",buflipbot," lipbufthick",lipbufthick + write (iout,*) "xi yi zi",xi,yi,zi +#endif if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then C the energy transfer exist if (zi.lt.buflipbot) then @@ -13275,5 +13394,8 @@ C lipbufthick is thickenes of lipid buffore sslipi=0.0d0 ssgradlipi=0.0 endif +#ifdef DEBUG + write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi +#endif return end diff --git a/source/unres/src-HCD-5D/energy_split-sep.F b/source/unres/src-HCD-5D/energy_split-sep.F index 11ea406..34a1bd1 100644 --- a/source/unres/src-HCD-5D/energy_split-sep.F +++ b/source/unres/src-HCD-5D/energy_split-sep.F @@ -118,10 +118,10 @@ C FG slaves receive the WEIGHTS array wliptran=weights(22) wtube=weights(25) wsaxs=weights(26) - wdfa_dist=weights_(28) - wdfa_tor=weights_(29) - wdfa_nei=weights_(30) - wdfa_beta=weights_(31) + wdfa_dist=weights(28) + wdfa_tor=weights(29) + wdfa_nei=weights(30) + wdfa_beta=weights(31) endif call MPI_Bcast(dc(1,1),6*nres,MPI_DOUBLE_PRECISION, & king,FG_COMM,IERR) diff --git a/source/unres/src-HCD-5D/gradient_p.F b/source/unres/src-HCD-5D/gradient_p.F index adafa53..67275ed 100644 --- a/source/unres/src-HCD-5D/gradient_p.F +++ b/source/unres/src-HCD-5D/gradient_p.F @@ -278,17 +278,20 @@ cd write(iout,*) 'calling int_to_cart' #ifdef DEBUG write (iout,*) "gcart, gxcart, gloc before int_to_cart" #endif - do i=1,nct + do i=0,nct do j=1,3 gcart(j,i)=gradc(j,i,icg) gxcart(j,i)=gradx(j,i,icg) enddo #ifdef DEBUG - if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + if (i.eq.0) then + write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3), + & (gxcart(j,i),j=1,3) + else if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg), & gloc(ialph(i,1),icg),gloc(ialph(i,1)+nside,icg) - else + else write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg) endif diff --git a/source/unres/src-HCD-5D/parmread.F b/source/unres/src-HCD-5D/parmread.F index 4da2913..7550fd5 100644 --- a/source/unres/src-HCD-5D/parmread.F +++ b/source/unres/src-HCD-5D/parmread.F @@ -69,6 +69,7 @@ C setenv LATEX YES C call getenv_loc("PRINT_PARM",lancuch) lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y") + & .and. (me.eq.king.or..not.out1file) .and. fg_rank.eq.0 call getenv_loc("LATEX",lancuch) LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y") C @@ -121,6 +122,7 @@ c & vbldsc0(j,i),aksc(j,i),abond0(j,i) enddo enddo + call flush(iout) endif C reading lipid parameters if (lprint) then @@ -132,6 +134,13 @@ C reading lipid parameters read(iliptranpar,*) liptranene(i) enddo close(iliptranpar) + if (lprint) then + write (iout,'(/a)') "Water-lipid transfer parameters" + write (iout,'(a3,3x,f10.5)') 'p',pepliptran + do i=1,ntyp + write (iout,'(a3,3x,f10.5)') restyp(i),liptranene(i) + enddo + endif #ifdef CRYST_THETA C C Read the parameters of the probability distribution/energy expression @@ -1452,7 +1461,7 @@ cc maxinter is maximum interaction sites #endif if (lprint) then - write (iout,'(/a/)') 'Torsional constants:' + write (iout,'(/a/)') 'SCCor torsional constants:' do l=1,maxinter do i=1,nsccortyp do j=1,nsccortyp @@ -1469,6 +1478,7 @@ cc maxinter is maximum interaction sites enddo enddo enddo + call flush(iout) endif C @@ -1951,6 +1961,7 @@ c call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) + call reada(weightcard,'LIPSCALE',lipscale,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) diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F index 8346c4a..41fe7f6 100644 --- a/source/unres/src-HCD-5D/readpdb-mult.F +++ b/source/unres/src-HCD-5D/readpdb-mult.F @@ -342,6 +342,7 @@ c write(iout,*)"before int_from_cart nres",nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo + dc(:,0)=c(:,1) do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) diff --git a/source/unres/src-HCD-5D/readpdb.F b/source/unres/src-HCD-5D/readpdb.F deleted file mode 100644 index c56b1df..0000000 --- a/source/unres/src-HCD-5D/readpdb.F +++ /dev/null @@ -1,631 +0,0 @@ - subroutine readpdb -C Read the PDB file and convert the peptide geometry into virtual-chain -C geometry. - implicit none - include 'DIMENSIONS' - include 'COMMON.LOCAL' - include 'COMMON.VAR' - include 'COMMON.CHAIN' - include 'COMMON.INTERACT' - include 'COMMON.IOUNITS' - include 'COMMON.GEO' - include 'COMMON.NAMES' - include 'COMMON.CONTROL' - include 'COMMON.FRAG' - include 'COMMON.SETUP' - include 'COMMON.SBRIDGE' - character*3 seq,atom,res - character*80 card - double precision sccor(3,50) - double precision e1(3),e2(3),e3(3) - integer rescode,iterter(maxres),cou - logical fail - integer i,j,iii,ires,ires_old,ishift,ibeg - double precision dcj - bfac=0.0d0 - do i=1,maxres - iterter(i)=0 - enddo - ibeg=1 - lsecondary=.false. - nhfrag=0 - nbfrag=0 - ires=0 - do - read (ipdbin,'(a80)',end=10) card - if (card(:5).eq.'HELIX') then - nhfrag=nhfrag+1 - lsecondary=.true. - read(card(22:25),*) hfrag(1,nhfrag) - read(card(34:37),*) hfrag(2,nhfrag) - endif - if (card(:5).eq.'SHEET') then - nbfrag=nbfrag+1 - lsecondary=.true. - read(card(24:26),*) bfrag(1,nbfrag) - read(card(35:37),*) bfrag(2,nbfrag) -crc---------------------------------------- -crc to be corrected !!! - bfrag(3,nbfrag)=bfrag(1,nbfrag) - bfrag(4,nbfrag)=bfrag(2,nbfrag) -crc---------------------------------------- - endif - if (card(:3).eq.'END') then - goto 10 - else if (card(:3).eq.'TER') then -C End current chain - ires_old=ires+2 - itype(ires_old-1)=ntyp1 - iterter(ires_old-1)=1 - itype(ires_old)=ntyp1 - iterter(ires_old)=1 - ibeg=2 - write (iout,*) "Chain ended",ires,ishift,ires_old - if (unres_pdb) then - do j=1,3 - dc(j,ires)=sccor(j,iii) - enddo - else - call sccenter(ires,iii,sccor) - endif - endif -C Fish out the ATOM cards. - if (index(card(1:4),'ATOM').gt.0) then - read (card(14:16),'(a3)') atom - if (atom.eq.'CA' .or. atom.eq.'CH3') then -C Calculate the CM of the preceding residue. - if (ibeg.eq.0) then - if (unres_pdb) then - do j=1,3 - dc(j,ires+nres)=sccor(j,iii) - enddo - else - call sccenter(ires,iii,sccor) - endif - endif -C Start new residue. -c write (iout,'(a80)') card - read (card(23:26),*) ires - read (card(18:20),'(a3)') res - if (ibeg.eq.1) then - ishift=ires-1 - if (res.ne.'GLY' .and. res.ne. 'ACE') then - ishift=ishift-1 - itype(1)=ntyp1 - endif -c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift - ibeg=0 - else if (ibeg.eq.2) then -c Start a new chain - ishift=-ires_old+ires-1 -c write (iout,*) "New chain started",ires,ishift - ibeg=0 - endif - ires=ires-ishift -c write (2,*) "ires",ires," ishift",ishift - if (res.eq.'ACE') then - itype(ires)=10 - else - itype(ires)=rescode(ires,res,0) - endif - read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) - read(card(61:66),*) bfac(ires) -c if(me.eq.king.or..not.out1file) -c & write (iout,'(2i3,2x,a,3f8.3)') -c & ires,itype(ires),res,(c(j,ires),j=1,3) - iii=1 - do j=1,3 - sccor(j,iii)=c(j,ires) - enddo - else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and. - & atom.ne.'N ' .and. atom.ne.'C ') then - iii=iii+1 - read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) - endif - endif - enddo - 10 if(me.eq.king.or..not.out1file) - & write (iout,'(a,i5)') ' Nres: ',ires -C Calculate dummy residue coordinates inside the "chain" of a multichain -C system - nres=ires - do i=2,nres-1 -c write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i) - if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then - if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then -C 16/01/2014 by Adasko: Adding to dummy atoms in the chain -C first is connected prevous chain (itype(i+1).eq.ntyp1)=true -C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the last dummy residue - print *,i,'tu dochodze' - call refsys(i-3,i-2,i-1,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif !fail - print *,i,'a tu?' - do j=1,3 - c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) - enddo - else !unres_pdb - do j=1,3 - dcj=(c(j,i-2)-c(j,i-3))/2.0 - if (dcj.eq.0) dcj=1.23591524223 - c(j,i)=c(j,i-1)+dcj - c(j,nres+i)=c(j,i) - enddo - endif !unres_pdb - else !itype(i+1).eq.ntyp1 - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the first dummy residue - call refsys(i+1,i+2,i+3,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,i)=c(j,i+1)-1.9d0*e2(j) - enddo - else !unres_pdb - do j=1,3 - dcj=(c(j,i+3)-c(j,i+2))/2.0 - if (dcj.eq.0) dcj=1.23591524223 - c(j,i)=c(j,i+1)-dcj - c(j,nres+i)=c(j,i) - enddo - endif !unres_pdb - endif !itype(i+1).eq.ntyp1 - endif !itype.eq.ntyp1 - enddo - write (iout,*) "After loop in readpbd" -C Calculate the CM of the last side chain. - if (unres_pdb) then - do j=1,3 - dc(j,ires)=sccor(j,iii) - enddo - else - call sccenter(ires,iii,sccor) - endif - nsup=nres - nstart_sup=1 - if (itype(nres).ne.10) then - nres=nres+1 - itype(nres)=ntyp1 - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the last dummy residue - call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) - enddo - else - do j=1,3 - dcj=(c(j,nres-2)-c(j,nres-3))/2.0 - if (dcj.eq.0) dcj=1.23591524223 - c(j,nres)=c(j,nres-1)+dcj - c(j,2*nres)=c(j,nres) - enddo - endif - endif - do i=2,nres-1 - do j=1,3 - c(j,i+nres)=dc(j,i) - enddo - enddo - do j=1,3 - c(j,nres+1)=c(j,1) - c(j,2*nres)=c(j,nres) - enddo - if (itype(1).eq.ntyp1) then - nsup=nsup-1 - nstart_sup=2 - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the first dummy residue - call refsys(2,3,4,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) - enddo - else - do j=1,3 - dcj=(c(j,4)-c(j,3))/2.0 - c(j,1)=c(j,2)-dcj - c(j,nres+1)=c(j,1) - enddo - endif - endif -C Calculate internal coordinates. - if(me.eq.king.or..not.out1file)then - do ires=1,nres - write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') - & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3), - & (c(j,nres+ires),j=1,3) - enddo - endif - call flush(iout) -c write(iout,*)"before int_from_cart nres",nres - call int_from_cart(.true.,.false.) - do i=1,nres - thetaref(i)=theta(i) - phiref(i)=phi(i) - enddo - 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 -c write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3), -c & vbld_inv(i+1) - enddo - do i=2,nres-1 - 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 -c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), -c & vbld_inv(i+nres) - enddo - call sc_loc_geom(.false.) - call int_from_cart1(.false.) -c call chainbuild -C Copy the coordinates to reference coordinates - do i=1,nres - do j=1,3 - cref(j,i)=c(j,i) - cref(j,i+nres)=c(j,i+nres) - enddo - enddo - 100 format (//' alpha-carbon coordinates ', - & ' centroid coordinates'/ - 1 ' ', 6X,'X',11X,'Y',11X,'Z', - & 10X,'X',11X,'Y',11X,'Z') - 110 format (a,'(',i3,')',6f12.5) -cc enddiag - do j=1,nbfrag - do i=1,4 - bfrag(i,j)=bfrag(i,j)-ishift - enddo - enddo - - do j=1,nhfrag - do i=1,2 - hfrag(i,j)=hfrag(i,j)-ishift - enddo - enddo - return - end -c--------------------------------------------------------------------------- - subroutine readpdb_template(k) -C Read the PDB file for read_constr_homology with read2sigma -C and convert the peptide geometry into virtual-chain geometry. - implicit none - include 'DIMENSIONS' - include 'COMMON.LOCAL' - include 'COMMON.VAR' - include 'COMMON.CHAIN' - include 'COMMON.INTERACT' - include 'COMMON.IOUNITS' - include 'COMMON.GEO' - include 'COMMON.NAMES' - include 'COMMON.CONTROL' - include 'COMMON.FRAG' - include 'COMMON.SETUP' - integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, - & ishift_pdb,ires_ca - logical lprn /.false./,fail - double precision e1(3),e2(3),e3(3) - double precision dcj,efree_temp - character*3 seq,res - character*5 atom - character*80 card - double precision sccor(3,20) - integer rescode,iterter(maxres) - do i=1,maxres - iterter(i)=0 - enddo - ibeg=1 - ishift1=0 - ishift=0 -c write (2,*) "UNRES_PDB",unres_pdb - ires=0 - ires_old=0 - iii=0 - lsecondary=.false. - nhfrag=0 - nbfrag=0 - do - read (ipdbin,'(a80)',end=10) card - if (card(:3).eq.'END') then - goto 10 - else if (card(:3).eq.'TER') then -C End current chain - ires_old=ires+2 - itype(ires_old-1)=ntyp1 - iterter(ires_old-1)=1 - itype(ires_old)=ntyp1 - iterter(ires_old)=1 - ibeg=2 -c write (iout,*) "Chain ended",ires,ishift,ires_old - if (unres_pdb) then - do j=1,3 - dc(j,ires)=sccor(j,iii) - enddo - else - call sccenter(ires,iii,sccor) - endif - endif -C Fish out the ATOM cards. - if (index(card(1:4),'ATOM').gt.0) then - read (card(12:16),*) atom -c write (iout,*) "! ",atom," !",ires -c if (atom.eq.'CA' .or. atom.eq.'CH3') then - read (card(23:26),*) ires - read (card(18:20),'(a3)') res -c write (iout,*) "ires",ires,ires-ishift+ishift1, -c & " ires_old",ires_old -c write (iout,*) "ishift",ishift," ishift1",ishift1 -c write (iout,*) "IRES",ires-ishift+ishift1,ires_old - if (ires-ishift+ishift1.ne.ires_old) then -C Calculate the CM of the preceding residue. - if (ibeg.eq.0) then - if (unres_pdb) then - do j=1,3 - dc(j,ires)=sccor(j,iii) - enddo - else - call sccenter(ires_old,iii,sccor) - endif - iii=0 - endif -C Start new residue. - if (res.eq.'Cl-' .or. res.eq.'Na+') then - ires=ires_old - cycle - else if (ibeg.eq.1) then -c write (iout,*) "BEG ires",ires - ishift=ires-1 - if (res.ne.'GLY' .and. res.ne. 'ACE') then - ishift=ishift-1 - itype(1)=ntyp1 - endif - ires=ires-ishift+ishift1 - ires_old=ires -c write (iout,*) "ishift",ishift," ires",ires, -c & " ires_old",ires_old -c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift - ibeg=0 - else if (ibeg.eq.2) then -c Start a new chain - ishift=-ires_old+ires-1 - ires=ires_old+1 -c write (iout,*) "New chain started",ires,ishift - ibeg=0 - else - ishift=ishift-(ires-ishift+ishift1-ires_old-1) - ires=ires-ishift+ishift1 - ires_old=ires - endif - if (res.eq.'ACE' .or. res.eq.'NHE') then - itype(ires)=10 - else - itype(ires)=rescode(ires,res,0) - endif - else - ires=ires-ishift+ishift1 - endif -c write (iout,*) "ires_old",ires_old," ires",ires -c if (card(27:27).eq."A" .or. card(27:27).eq."B") then -c ishift1=ishift1+1 -c endif -c write (2,*) "ires",ires," res ",res," ity",ity - if (atom.eq.'CA' .or. atom.eq.'CH3' .or. - & res.eq.'NHE'.and.atom(:2).eq.'HN') then - read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) -c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3) -#ifdef DEBUG - write (iout,'(2i3,2x,a,3f8.3)') - & ires,itype(ires),res,(c(j,ires),j=1,3) -#endif - iii=iii+1 - do j=1,3 - sccor(j,iii)=c(j,ires) - enddo - if (ishift.ne.0) then - ires_ca=ires+ishift-ishift1 - else - ires_ca=ires - endif -c write (*,*) card(23:27),ires,itype(ires) - else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. - & atom.ne.'N' .and. atom.ne.'C' .and. - & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. - & atom.ne.'OXT' .and. atom(:2).ne.'3H') then -c write (iout,*) "sidechain ",atom - iii=iii+1 - read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) - endif - endif - enddo - 10 if(me.eq.king.or..not.out1file) - & write (iout,'(a,i5)') ' Nres: ',ires -C Calculate dummy residue coordinates inside the "chain" of a multichain -C system - nres=ires - do i=2,nres-1 -c write (iout,*) i,itype(i),itype(i+1) - if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then - if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then -C 16/01/2014 by Adasko: Adding to dummy atoms in the chain -C first is connected prevous chain (itype(i+1).eq.ntyp1)=true -C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the last dummy residue - call refsys(i-3,i-2,i-1,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif !fail - do j=1,3 - c(j,i)=c(j,i-1)-1.9d0*e2(j) - enddo - else !unres_pdb - do j=1,3 - dcj=(c(j,i-2)-c(j,i-3))/2.0 - if (dcj.eq.0) dcj=1.23591524223 - c(j,i)=c(j,i-1)+dcj - c(j,nres+i)=c(j,i) - enddo - endif !unres_pdb - else !itype(i+1).eq.ntyp1 - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the first dummy residue - call refsys(i+1,i+2,i+3,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,i)=c(j,i+1)-1.9d0*e2(j) - enddo - else !unres_pdb - do j=1,3 - dcj=(c(j,i+3)-c(j,i+2))/2.0 - if (dcj.eq.0) dcj=1.23591524223 - c(j,i)=c(j,i+1)-dcj - c(j,nres+i)=c(j,i) - enddo - endif !unres_pdb - endif !itype(i+1).eq.ntyp1 - endif !itype.eq.ntyp1 - enddo -C Calculate the CM of the last side chain. - if (unres_pdb) then - do j=1,3 - dc(j,ires)=sccor(j,iii) - enddo - else - call sccenter(ires,iii,sccor) - endif - nsup=nres - nstart_sup=1 - if (itype(nres).ne.10) then - nres=nres+1 - itype(nres)=ntyp1 - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the last dummy residue - call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,nres)=c(j,nres-1)-1.9d0*e2(j) - enddo - else - do j=1,3 - dcj=(c(j,nres-2)-c(j,nres-3))/2.0 - if (dcj.eq.0) dcj=1.23591524223 - c(j,nres)=c(j,nres-1)+dcj - c(j,2*nres)=c(j,nres) - enddo - endif - endif - do i=2,nres-1 - do j=1,3 - c(j,i+nres)=dc(j,i) - enddo - enddo - do j=1,3 - c(j,nres+1)=c(j,1) - c(j,2*nres)=c(j,nres) - enddo - if (itype(1).eq.ntyp1) then - nsup=nsup-1 - nstart_sup=2 - if (unres_pdb) then -C 2/15/2013 by Adam: corrected insertion of the first dummy residue - call refsys(2,3,4,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,1)=c(j,2)-1.9d0*e2(j) - enddo - else - do j=1,3 - dcj=(c(j,4)-c(j,3))/2.0 - c(j,1)=c(j,2)-dcj - c(j,nres+1)=c(j,1) - enddo - endif - endif -C Copy the coordinates to reference coordinates -c do i=1,2*nres -c do j=1,3 -c cref(j,i)=c(j,i) -c enddo -c enddo -C Calculate internal coordinates. - if (out_template_coord) then - write (iout,'(/a)') - & "Cartesian coordinates of the reference structure" - write (iout,'(a,3(3x,a5),5x,3(3x,a5))') - & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" - do ires=1,nres - write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)') - & restyp(itype(ires)),ires,(c(j,ires),j=1,3), - & (c(j,ires+nres),j=1,3) - enddo - endif -C Calculate internal coordinates. - call int_from_cart(.true.,.true.) - call sc_loc_geom(.false.) - do i=1,nres - thetaref(i)=theta(i) - phiref(i)=phi(i) - enddo - 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=2,nres-1 - 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 -c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), -c & vbld_inv(i+nres) - enddo - do i=1,nres - do j=1,3 - cref(j,i)=c(j,i) - cref(j,i+nres)=c(j,i+nres) - enddo - enddo - do i=1,2*nres - do j=1,3 - chomo(j,i,k)=c(j,i) - enddo - enddo - - return - end - diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index f120ec7..da28aa3 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -1215,8 +1215,6 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb C endif -c write (iout,*) "iranconf",iranconf," extconf",extconf, -c & " start_from_models",start_from_model 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 @@ -1230,11 +1228,16 @@ C initial geometry. 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) + if (nnt.gt.1) c(:,nres+1)=c(:,1) + if (nct.lt.nres) c(:,2*nres)=c(:,nres) c write (iout,*) "Exit READ_CART" c write (iout,'(8f10.5)') c & ((c(l,k),l=1,3),k=1,nres), c & ((c(l,k+nres),l=1,3),k=nnt,nct) call cartprint + do j=1,3 + dc(j,0)=c(j,1) + enddo do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) diff --git a/source/unres/src-HCD-5D/unres.F b/source/unres/src-HCD-5D/unres.F index be40e75..978fa59 100644 --- a/source/unres/src-HCD-5D/unres.F +++ b/source/unres/src-HCD-5D/unres.F @@ -59,7 +59,8 @@ c call memmon_print_usage() C Read force field parameters and job setup data call readrtns C -c write (iout,*) "After readrtns" + write (iout,*) "After readrtns" + call flush(iout) call cartprint call intout if (me.eq.king .or. .not. out1file) then @@ -86,9 +87,14 @@ C Fine-grain slaves just do energy and gradient components. call ergastulum ! slave workhouse in Latin else #endif + if (indpdb.eq.0 .and. .not.read_cart) call chainbuild + if (indpdb.ne.0 .or. read_cart) then + dc(1,0)=c(1,1) + dc(2,0)=c(2,1) + dc(3,0)=c(3,1) + endif if (modecalc.eq.0) then write (iout,*) "Calling exec_eeval_or_minim" - call cartprint call exec_eeval_or_minim else if (modecalc.eq.1) then call exec_regularize @@ -223,8 +229,8 @@ c include 'COMMON.CONTACTS' common /lbfgstat/ status,niter,nfun #endif integer ilen - if (indpdb.eq.0) call chainbuild - if (indpdb.ne.0) then + if (indpdb.eq.0 .and. .not.read_cart) call chainbuild + if (indpdb.ne.0 .or. read_cart) then dc(1,0)=c(1,1) dc(2,0)=c(2,1) dc(3,0)=c(3,1) @@ -331,16 +337,7 @@ c print *,'Calling MINIMIZE.' call etotal(energy(0)) etot = energy(0) call enerprint(energy(0)) - call intout - if (out_int) call briefout(0,etot) - if (out_cart) then - cartname=prefix(:ilen(prefix))//'.x' - potE=etot - call cartoutx(0.0d0) - endif - if (outpdb) call pdbout(etot,titel(:50),ipdb) - if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) #ifdef LBFGS write (iout,'(a,a9)') 'LBFGS return code:',status write (iout,'(a,i20)') '# of energy evaluations:',nfun+1 @@ -350,10 +347,13 @@ c print *,'Calling MINIMIZE.' write (iout,'(a,i20)') '# of energy evaluations:',nfun+1 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals #endif - else - print *,'refstr=',refstr - if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) - call briefout(0,etot) + endif + if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + if (out_int) call briefout(0,etot) + if (out_cart) then + cartname=prefix(:ilen(prefix))//'.x' + potE=etot + call cartoutx(0.0d0) endif if (outpdb) call pdbout(etot,titel(:50),ipdb) if (outmol2) call mol2out(etot,titel(:32)) diff --git a/source/wham/src-HCD/Makefile b/source/wham/src-HCD/Makefile index ee054bf..8453cdd 120000 --- a/source/wham/src-HCD/Makefile +++ b/source/wham/src-HCD/Makefile @@ -1 +1 @@ -Makefile_MPICH_ifort-okeanos \ No newline at end of file +Makefile_MPICH_ifort \ No newline at end of file diff --git a/source/wham/src-HCD/Makefile_MPICH_ifort b/source/wham/src-HCD/Makefile_MPICH_ifort index 9a83c35..a04725d 100644 --- a/source/wham/src-HCD/Makefile_MPICH_ifort +++ b/source/wham/src-HCD/Makefile_MPICH_ifort @@ -1,10 +1,13 @@ -INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh +INSTALL_DIR = /users/software/mpich2-1.4.1p1_intel +OPT = -mcmodel=medium -shared-intel -O3 -dynamic BIN = ../../../bin/wham -FC= ifort -OPT = -mcmodel=medium -O3 -ip -w -#OPT = -mcmodel=medium -g -CB +FC= ${INSTALL_DIR}/bin/mpif90 +OPT = -mcmodel=medium -O3 -ip -w +#OPT = -O3 -intel-static -mcmodel=medium +#OPT = -O3 -ip -w +#OPT = -g -CB -mcmodel=medium -shared-intel -dynamic FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include -LIBS = -L$(INSTALL_DIR)/lib -lmpich -lpmpich xdrf/libxdrf.a +LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a .f.o: ${FC} ${FFLAGS} $*.f @@ -24,20 +27,24 @@ objects = \ molread_zs.o \ openunits.o \ readrtns.o \ + read_constr_homology.o \ arcos.o \ - cartder.o \ cartprint.o \ chainbuild.o \ geomout.o \ icant.o \ intcor.o \ int_from_cart.o \ + refsys.o \ make_ensemble1.o \ matmult.o \ misc.o \ mygetenv.o \ parmread.o \ permut.o \ + seq2chains.o \ + chain_symmetry.o \ + iperm.o \ pinorm.o \ printmat.o \ proc_proc.o \ @@ -47,7 +54,10 @@ objects = \ store_parm.o \ timing.o \ wham_calc1.o \ - ssMD.o + PMFprocess.o \ + ssMD.o \ + boxshift.o \ + oligomer.o objects_compar = \ readrtns_compar.o \ @@ -57,43 +67,93 @@ objects_compar = \ rmscalc.o secondary.o proc_cont.o define_pairs.o mysort.o all: no_option - @echo "Specify force field: GAB, 4P or E0LL2Y" + @echo "Specify force field: GAB, 4P, E0LL2Y or NEWCORR" no_option: GAB: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DAMD64 \ - -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DWHAM + -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DFOURBODY -DWHAM GAB: ${objects} ${objects_compar} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c ./compinfo ${FC} -c ${FFLAGS} cinfo.f $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \ - ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH_GAB.exe + ${LIBS} -o ${BIN}/wham_ifort_MPICH-GAB-HCD.exe + +GAB_DFA: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DAMD64 \ + -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DFOURBODY -DWHAM -DDFA +GAB_DFA: ${objects} ${objects_compar} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo + ${FC} -c ${FFLAGS} cinfo.f + $(FC) ${OPT} ${objects} ${objects_compar} dfa.o cinfo.o \ + ${LIBS} -o ${BIN}/wham_ifort_MPICH-GAB-HCD-DFA.exe 4P: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPGI -DISNAN -DAMD64 \ - -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DWHAM + -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DFOURBODY -DWHAM 4P: ${objects} ${objects_compar} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c ./compinfo ${FC} -c ${FFLAGS} cinfo.f $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \ - ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH_4P.exe + ${LIBS} -o ${BIN}/wham_ifort_MPICH-4P-HCD.exe + +4P_DFA: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPGI -DISNAN -DAMD64 \ + -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DFOURBODY -DWHAM -DDFA +4P_DFA: ${objects} ${objects_compar} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo + ${FC} -c ${FFLAGS} cinfo.f + $(FC) ${OPT} ${objects} ${objects_compar} dfa.o cinfo.o \ + ${LIBS} -o ${BIN}/wham_ifort_MPICH-4P-HCD-DFA.exe -E0LL2Y: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DAMD64 -DWHAM +E0LL2Y: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DFOURBODY -DAMD64 -DWHAM E0LL2Y: ${objects} ${objects_compar} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c ./compinfo ${FC} -c ${FFLAGS} cinfo.f $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \ - ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH_E0LL2Y.exe + ${LIBS} -o ${BIN}/wham_ifort_MPICH-E0LL2Y-HCD.exe -NEWCORR: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DNEWCORR -DPGI -DISNAN -DAMD64 -DWHAM +E0LL2Y_DFA: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DFOURBODY -DAMD64 -DWHAM -DDFA +E0LL2Y_DFA: ${objects} ${objects_compar} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo + ${FC} -c ${FFLAGS} cinfo.f + $(FC) ${OPT} ${objects} ${objects_compar} dfa.o cinfo.o \ + ${LIBS} -o ${BIN}/wham_ifort_MPICH-E0LL2Y-HCD-DFA.exe + +NEWCORR: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DNEWCORR -DCORRCD -DPGI -DISNAN -DAMD64 -DWHAM NEWCORR: ${objects} ${objects_compar} xdrf/libxdrf.a - cc -o compinfo compinfo.c + gcc -o compinfo compinfo.c ./compinfo ${FC} -c ${FFLAGS} cinfo.f $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \ - ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH_NEWCORR.exe + ${LIBS} -o ${BIN}/wham_ifort_MPICH-SC-HCD.exe + +NEWCORR5D: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DNEWCORR -DCORRCD -DPGI -DISNAN -DAMD64 -DWHAM -DFIVEDIAG +NEWCORR5D: ${objects} ${objects_compar} xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo + ${FC} -c ${FFLAGS} cinfo.f + $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \ + ${LIBS} -o ${BIN}/wham_ifort_MPICH-SC-HCD5.exe + +NEWCORR_DFA: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DNEWCORR -DCORRCD -DDFA -DPGI -DISNAN -DAMD64 -DWHAM -DDFA +NEWCORR_DFA: ${objects} ${objects_compar} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo + ${FC} -c ${FFLAGS} cinfo.f + $(FC) ${OPT} ${objects} ${objects_compar} dfa.o cinfo.o \ + ${LIBS} -o ${BIN}/wham_ifort_MPICH-SC-HCD-DFA-D.exe + +NEWCORR5D_DFA: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DNEWCORR -DCORRCD -DDFA -DPGI -DISNAN -DAMD64 -DWHAM -DFIVEDIAG -DDFA +NEWCORR5D_DFA: ${objects} ${objects_compar} dfa.o xdrf/libxdrf.a + gcc -o compinfo compinfo.c + ./compinfo + ${FC} -c ${FFLAGS} cinfo.f + $(FC) ${OPT} ${objects} ${objects_compar} dfa.o cinfo.o \ + ${LIBS} -o ${BIN}/wham_ifort_MPICH-SC-HCD5-DFA.exe xdrf/libxdrf.a: cd xdrf && make diff --git a/source/wham/src-HCD/energy_p_new.F b/source/wham/src-HCD/energy_p_new.F index ce7a6a7..e72d558 100644 --- a/source/wham/src-HCD/energy_p_new.F +++ b/source/wham/src-HCD/energy_p_new.F @@ -1216,6 +1216,9 @@ C returning jth atom to box & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 +c write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) +c if (aa.ne.aa_aq(itypi,itypj)) write(iout,'(2e15.5)') +c &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1278,8 +1281,8 @@ c#define DEBUG #endif c#undef DEBUG c endif - if (energy_dec) write (iout,'(a,2i5,3f10.5)') - & 'r sss evdw',i,j,1.0d0/rij,sss,evdwij + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij if (calc_grad) then C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1288,6 +1291,12 @@ C Calculate gradient components. fac=rij*fac fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -1457,6 +1466,12 @@ C Calculate gradient components. fac=rij*fac-2*expon*rrij*e_augm fac=fac+(evdwij+e_augm)*sssgrad/sss*rij C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -2110,6 +2125,8 @@ C common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -2217,6 +2234,7 @@ c end if ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -2247,6 +2265,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2287,6 +2306,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2358,6 +2378,9 @@ C------------------------------------------------------------------------------- common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij, + & faclipij2 + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -2393,6 +2416,9 @@ C zj=c(3,j)+0.5D0*dzj-zmedi yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0 + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -2432,25 +2458,25 @@ C fac_shield(j)=0.6 el1=el1*fac_shield(i)**2*fac_shield(j)**2 el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*faclipij2 else fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*faclipij2 endif - evdw1=evdw1+evdwij*sss + evdw1=evdw1+evdwij*sss*faclipij2 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') - &'evdw1',i,j,evdwij - &,iteli,itelj,aaa,evdw1,sss - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - &fac_shield(i),fac_shield(j) + write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') + &' evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss + write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij, + & faclipij2 endif C @@ -2468,9 +2494,10 @@ C * Radial derivatives. First process both termini of the fragment (i,j) * if (calc_grad) then - ggg(1)=facel*xj - ggg(2)=facel*yj - ggg(3)=facel*zj + aux=(facel*sss+rmij*sssgrad*eesij)*faclipij2 + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j @@ -2556,6 +2583,11 @@ C gelc_long(k,i-1)=gelc_long(k,i-1) C & +grad_shield(k,i)*eesij/fac_shield(i) C gelc_long(k,j-1)=gelc_long(k,j-1) C & +grad_shield(k,j)*eesij/fac_shield(j) + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2*sss + + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2*sss enddo C print *,"bafter", gelc_long(1,i), gelc_long(1,j) @@ -2568,7 +2600,7 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo if (sss.gt.0.0) then - facvdw=facvdw+sssgrad*rmij*evdwij + facvdw=(facvdw+sssgrad*rmij*evdwij)*faclipij2 ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj @@ -2587,6 +2619,11 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +!C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 * * Loop over residues i+1 thru j-1. * @@ -2598,7 +2635,7 @@ cgrad enddo endif ! calc_grad #else C MARYSIA - facvdw=(ev1+evdwij) + facvdw=(ev1+evdwij)*faclipij2 facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel)*sss @@ -2642,6 +2679,10 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 endif ! calc_grad #endif * @@ -2661,7 +2702,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* - & fac_shield(i)**2*fac_shield(j)**2 + & fac_shield(i)**2*fac_shield(j)**2*sss*faclipij2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -2682,11 +2723,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -2918,7 +2959,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij c if (eel_loc_ij.ne.0) @@ -2982,7 +3023,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -2998,7 +3039,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -3011,7 +3052,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -3023,7 +3064,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -3032,12 +3073,12 @@ C Partial derivatives in virtual-bond dihedral angles gamma & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) aux=eel_loc_ij/sss*sssgrad*rmij ggg(1)=aux*xj @@ -3046,13 +3087,18 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/faclipij + + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -3062,19 +3108,19 @@ C Remaining derivatives of eello do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij enddo endif ! calc_grad @@ -3342,6 +3388,8 @@ C Third- and fourth-order contributions from turns common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij j=i+2 c write (iout,*) "eturn3",i,j,j1,j2 a_temp(1,1)=a22 @@ -3379,7 +3427,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) if (energy_dec) write (iout,'(6heturn3,2i5,0pf7.3)') i,i+2, @@ -3389,10 +3437,10 @@ C#ifdef NEWCORR C Derivatives in theta gloc(nphi+i,icg)=gloc(nphi+i,icg) & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C#endif C Derivatives in shield mode @@ -3447,14 +3495,14 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Cartesian derivatives do l=1,3 c ghalf1=0.5d0*agg(l,1) @@ -3468,7 +3516,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) @@ -3477,7 +3525,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -3485,7 +3533,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -3493,7 +3541,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif ! calc_grad @@ -3628,7 +3676,7 @@ C fac_shield(i)=0.6 C fac_shield(j)=0.4 endif eello_turn4=eello_turn4-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) @@ -3705,7 +3753,7 @@ C Derivatives in gamma(i) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -3714,7 +3762,7 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -3726,7 +3774,7 @@ C Derivatives in gamma(i+2) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij if (calc_grad) then C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) @@ -3747,7 +3795,7 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif C Remaining derivatives of this turn contribution @@ -3766,7 +3814,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -3781,7 +3829,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -3796,7 +3844,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -3812,7 +3860,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif ! calc_grad @@ -9492,6 +9540,8 @@ C--bufliptop--- here true lipid starts C lipid C--buflipbot--- lipid ends buffore starts C--bordlipbot--buffore ends +c call cartprint +c write (iout,*) "Eliptransfer peplipran",pepliptran eliptran=0.0 do i=1,nres C do i=1,1 @@ -9543,6 +9593,8 @@ CV do i=1,1 if (itype(i).eq.ntyp1) cycle positi=(mod(c(3,i+nres),boxzsize)) if (positi.le.0) positi=positi+boxzsize +c write(iout,*) "i",i," positi",positi,bordlipbot,buflipbot, +c & bordliptop C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop c for each residue check if it is in lipid or lipid water border area C respos=mod(c(3,i+nres),boxzsize) @@ -9553,8 +9605,11 @@ C the energy transfer exist if (positi.lt.buflipbot) then fracinbuf=1.0d0- & ((positi-bordlipbot)/lipbufthick) +c write (iout,*) "i",i,itype(i)," fracinbuf",fracinbuf +c write (iout,*) "i",i," liptranene",liptranene(itype(i)) C lipbufthick is thickenes of lipid buffore sslip=sscalelip(fracinbuf) +c write (iout,*) "sslip",sslip ssgradlip=-sscagradlip(fracinbuf)/lipbufthick eliptran=eliptran+sslip*liptranene(itype(i)) gliptranx(3,i)=gliptranx(3,i) @@ -9580,6 +9635,7 @@ C print *,"I am in true lipid" endif ! if in lipid or buffor C else C eliptran=elpitran+0.0 ! I am in water +c write (iout,*) "eliptran",eliptran enddo return end diff --git a/source/wham/src-HCD/include_unres/COMMON.INTERACT b/source/wham/src-HCD/include_unres/COMMON.INTERACT index 7d6b59f..c929acb 100644 --- a/source/wham/src-HCD/include_unres/COMMON.INTERACT +++ b/source/wham/src-HCD/include_unres/COMMON.INTERACT @@ -30,7 +30,9 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters. & distchainmax,nbondterm(ntyp) &,vbldpDUM C 01/29/15 Lipidic parameters - double precision pepliptran,liptranene - common /lipid/ pepliptran,liptranene(ntyp) + double precision pepliptran,liptranene,lipscale,tubetranene, + & tubetranenepep + common /lipid/ pepliptran,liptranene(ntyp),lipscale + common /tubepar/ tubetranene(ntyp), tubetranenepep diff --git a/source/wham/src-HCD/initialize_p.F b/source/wham/src-HCD/initialize_p.F index a2281e5..8217663 100644 --- a/source/wham/src-HCD/initialize_p.F +++ b/source/wham/src-HCD/initialize_p.F @@ -304,7 +304,7 @@ c------------------------------------------------------------------------- ! 15 16 17 18 19 20 21 & "WHPB ","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC", ! 22 23 24 25 26 27 28 - & "WLIPTRAN","WAFM","WTHETC","WSHIELD","WSAXS","WHOMO","WDFAD", + & "WLT ","WAFM","WTHETC","WSHIELD","WSAXS","WHOMO","WDFAD", ! 29 30 31 & "WDFAT","WDFAN","WDFAB"/ data ww0 / diff --git a/source/wham/src-HCD/oligomer.F b/source/wham/src-HCD/oligomer.F index c63ea47..44dddcf 100644 --- a/source/wham/src-HCD/oligomer.F +++ b/source/wham/src-HCD/oligomer.F @@ -40,7 +40,8 @@ c print *,i,c(:,i) sumodl_min=1.0d10 do i=0,3**nchain-1 ii=i - do ichain=1,nchain + lshift(1)=0 + do ichain=2,nchain lshift(ichain)=mod(ii,3)-1 ii=ii/3 enddo @@ -133,6 +134,7 @@ c coordinate system enddo #endif do ichain=1,nchain +c write (iout,*) "shift_coord",shift_coord(:,ichain) do i=chain_border1(1,ichain),chain_border1(2,ichain) do j=1,3 c(j,i)=c(j,i)+shift_coord(j,ichain) diff --git a/source/wham/src-HCD/parmread.F b/source/wham/src-HCD/parmread.F index b21acb2..4636b60 100644 --- a/source/wham/src-HCD/parmread.F +++ b/source/wham/src-HCD/parmread.F @@ -63,6 +63,8 @@ C Assign virtual-bond length call reada(controlcard,'SCALSCP',scalscp,1.0d0) call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0) call reada(controlcard,'DELT_CORR',delt_corr,0.5d0) + call reada(controlcard,'LIPSCALE',lipscale,1.0D0) +c write (iout,*) "lipscale",lipscale r0_corr=cutoff_corr-delt_corr write (iout,*) "iparm",iparm," myparm",myparm diff --git a/source/wham/src-HCD/proc_cont.f b/source/wham/src-HCD/proc_cont.f index 9269496..c70652a 100644 --- a/source/wham/src-HCD/proc_cont.f +++ b/source/wham/src-HCD/proc_cont.f @@ -15,10 +15,10 @@ include 'COMMON.GEO' write (iout,*) "proc_cont: nlevel",nlevel if (nlevel.lt.0) then - write (iout,*) "call define_fragments" +c write (iout,*) "call define_fragments" call define_fragments else - write (iout,*) "call secondary2" +c write (iout,*) "call secondary2" call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref, & isec_ref) endif -- 1.7.9.5