Adam's lipid and dfa corrections
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Sun, 7 Jun 2020 10:15:53 +0000 (12:15 +0200)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Sun, 7 Jun 2020 10:15:53 +0000 (12:15 +0200)
32 files changed:
source/cluster/wham/src-HCD/Makefile
source/cluster/wham/src-HCD/Makefile-MPICH-ifort
source/cluster/wham/src-HCD/Makefile-MPICH-ifort-prometheus
source/cluster/wham/src-HCD/Makefile-tryton
source/cluster/wham/src-HCD/energy_p_new.F
source/cluster/wham/src-HCD/include_unres/COMMON.INTERACT
source/cluster/wham/src-HCD/initialize.f [deleted file]
source/cluster/wham/src-HCD/readrtns.F
source/cluster/wham/src-HCD/sizesclu.dat
source/unres/src-HCD-5D/COMMON.DFA
source/unres/src-HCD-5D/COMMON.INTERACT
source/unres/src-HCD-5D/DIMENSIONS
source/unres/src-HCD-5D/MD_A-MTS.F
source/unres/src-HCD-5D/checkder_p.F
source/unres/src-HCD-5D/dfa.F
source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
source/unres/src-HCD-5D/energy_p_new_barrier.F
source/unres/src-HCD-5D/energy_split-sep.F
source/unres/src-HCD-5D/gradient_p.F
source/unres/src-HCD-5D/parmread.F
source/unres/src-HCD-5D/readpdb-mult.F
source/unres/src-HCD-5D/readpdb.F [deleted file]
source/unres/src-HCD-5D/readrtns_CSA.F
source/unres/src-HCD-5D/unres.F
source/wham/src-HCD/Makefile
source/wham/src-HCD/Makefile_MPICH_ifort
source/wham/src-HCD/energy_p_new.F
source/wham/src-HCD/include_unres/COMMON.INTERACT
source/wham/src-HCD/initialize_p.F
source/wham/src-HCD/oligomer.F
source/wham/src-HCD/parmread.F
source/wham/src-HCD/proc_cont.f

index 8aee570..693492e 120000 (symlink)
@@ -1 +1 @@
-Makefile-MPICH-ifort-okeanos
\ No newline at end of file
+Makefile-MPICH-ifort
\ No newline at end of file
index 79b8d0f..907ce62 100644 (file)
@@ -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
 
index 1492755..e1b8f32 100644 (file)
@@ -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
 
index e887bc9..ec6ae53 100644 (file)
@@ -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
index db2e043..4fa79c5 100644 (file)
@@ -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
index 1c0b8db..a02f7e4 100644 (file)
@@ -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 (file)
index 12ea156..0000000
+++ /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 
index 057f1ac..e9e576f 100644 (file)
@@ -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)
index cb8572c..7d0d666 100644 (file)
@@ -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.
index 67a1a0d..6759f8b 100644 (file)
@@ -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
index 8c4876d..14416ad 100644 (file)
@@ -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
 
index 137b45d..9803b23 100644 (file)
@@ -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)
index fcef69e..e504cbd 100644 (file)
@@ -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)
index 48eedda..e1448db 100644 (file)
@@ -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
index f69b81a..62e8892 100644 (file)
@@ -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
index 28ba1d1..c4e54bc 100644 (file)
@@ -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
index 2c701ca..6d5b25f 100644 (file)
@@ -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
       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
             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
       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
             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
       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
index 11ea406..34a1bd1 100644 (file)
@@ -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)
index adafa53..67275ed 100644 (file)
@@ -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
index 4da2913..7550fd5 100644 (file)
@@ -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)
index 8346c4a..41fe7f6 100644 (file)
@@ -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 (file)
index c56b1df..0000000
+++ /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
-      
index f120ec7..da28aa3 100644 (file)
@@ -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)
index be40e75..978fa59 100644 (file)
@@ -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))
index ee054bf..8453cdd 120000 (symlink)
@@ -1 +1 @@
-Makefile_MPICH_ifort-okeanos
\ No newline at end of file
+Makefile_MPICH_ifort
\ No newline at end of file
index 9a83c35..a04725d 100644 (file)
@@ -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
index ce7a6a7..e72d558 100644 (file)
@@ -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
index 7d6b59f..c929acb 100644 (file)
@@ -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
 
 
index a2281e5..8217663 100644 (file)
@@ -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 /
index c63ea47..44dddcf 100644 (file)
@@ -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)
index b21acb2..4636b60 100644 (file)
@@ -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
index 9269496..c70652a 100644 (file)
       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