-Makefile-MPICH-ifort-okeanos
\ No newline at end of file
+Makefile-MPICH-ifort
\ No newline at end of file
-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
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
-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
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
$(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
$(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
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
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include "DIMENSIONS.COMPAR"
+ include 'COMMON.CONTROL'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.LOCAL'
#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
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
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
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/
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)
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
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
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/
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)
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
* 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
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
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
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)
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
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.
& +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)
& +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)
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)
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
& 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
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 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
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
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,
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
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)
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)
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
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)
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
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)
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))
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))
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)
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
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)
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)
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)
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
& 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
+++ /dev/null
- 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
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)
* 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.
& 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
& 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
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)
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)
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)
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
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!!
integer ica1, ica2,ica3,ica4,ica5
integer ishell, inca, itmp,iitmp
double precision wtmp
+ logical lprn /.false./
C
C READ DISTANCE
C
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
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)
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
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)
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
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
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)
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
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)
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=e1+e2
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
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)
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
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)
& 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
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
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)
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
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)
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
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
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
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)
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)
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)
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)
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.
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)
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)
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)
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)
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
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)
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)
& 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
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
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
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)
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)
& 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
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
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.
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)
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
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.
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)
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
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)
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/
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)
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
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)
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/
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)
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
*
* 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
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
*
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
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.
*
*
* 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
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
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
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)
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.
& +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)
& +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)
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)
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)
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)
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)
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)
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
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,
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)
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)
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
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
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
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
#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
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
#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
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
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)
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)
#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
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)
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)
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
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)
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)
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
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)
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
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)
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)
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)
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
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
& +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)
& 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
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
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)
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
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/
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)
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
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
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/
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
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)
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,
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
*
* 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
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.
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
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.
*
cgrad enddo
#else
C MARYSIA
- facvdw=(ev1+evdwij)
+ facvdw=(ev1+evdwij)*faclipij2
facel=(el1+eesij)
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)*sss
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
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)
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
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
& +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)
& +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)
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)
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
& 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
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)
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
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
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,
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
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)
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)
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
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)
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-------------------------------------------------------------------------------
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)
& 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
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))
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))
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
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
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)
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)
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)
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-----------------------------------------------------------------------------
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
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)
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
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
sslipi=0.0d0
ssgradlipi=0.0
endif
+#ifdef DEBUG
+ write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi
+#endif
return
end
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)
#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
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
& vbldsc0(j,i),aksc(j,i),abond0(j,i)
enddo
enddo
+ call flush(iout)
endif
C reading lipid parameters
if (lprint) then
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
#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
enddo
enddo
enddo
+ call flush(iout)
endif
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)
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)
+++ /dev/null
- 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
-
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
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)
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
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
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)
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
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))
-Makefile_MPICH_ifort-okeanos
\ No newline at end of file
+Makefile_MPICH_ifort
\ No newline at end of file
-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
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 \
store_parm.o \
timing.o \
wham_calc1.o \
- ssMD.o
+ PMFprocess.o \
+ ssMD.o \
+ boxshift.o \
+ oligomer.o
objects_compar = \
readrtns_compar.o \
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
& +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)
#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
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
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
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/
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)
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
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
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/
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)
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
* 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
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)
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
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.
*
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
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
*
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)
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
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)
& +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)
& +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)
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)
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
& 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
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)
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
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
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,
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
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)
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)
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
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)
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
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)
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))
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))
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)
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
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)
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)
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)
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
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
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)
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)
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
& 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
! 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 /
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
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)
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
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
projects / unres.git / commitdiff