subroutine etotal(energia)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifndef ISNAN
external proc_proc
#ifdef MPI
include "mpif.h"
double precision weights_(n_ene)
+ double precision time00
+ integer ierror,ierr
#endif
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.VAR'
- include 'COMMON.MD'
+c include 'COMMON.MD'
+ include 'COMMON.QRESTR'
include 'COMMON.CONTROL'
include 'COMMON.TIME1'
include 'COMMON.SPLITELE'
include 'COMMON.TORCNSTR'
+ include 'COMMON.SAXS'
+ include 'COMMON.MD'
+ double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
+ & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
+ & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6,
+ & eliptran,Eafmforce,Etube,
+ & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
+ integer n_corr,n_corr1
#ifdef MPI
c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
c & " nfgtasks",nfgtasks
weights_(17)=wbond
weights_(18)=scal14
weights_(21)=wsccor
- weights_(22)=wtube
+ weights_(22)=wliptran
+ weights_(25)=wtube
weights_(26)=wsaxs
weights_(28)=wdfa_dist
weights_(29)=wdfa_tor
wbond=weights(17)
scal14=weights(18)
wsccor=weights(21)
- wtube=weights(22)
+ wliptran=weights(22)
+ wtube=weights(25)
wsaxs=weights(26)
wdfa_dist=weights_(28)
wdfa_tor=weights_(29)
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
c call chainbuild_cart
endif
-#ifndef DFA
- edfadis=0.0d0
- edfator=0.0d0
- edfanei=0.0d0
- edfabet=0.0d0
-#endif
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+ endif
+c write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate
+ if (mod(itime_mat,imatupdate).eq.0) then
+ call make_SCp_inter_list
+ call make_SCSC_inter_list
+ call make_pp_inter_list
+ call make_pp_vdw_inter_list
+ endif
c print *,'Processor',myrank,' calling etotal ipot=',ipot
c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
#else
#ifdef TIMING
time00=MPI_Wtime()
#endif
+
+#ifndef DFA
+ edfadis=0.0d0
+ edfator=0.0d0
+ edfanei=0.0d0
+ edfabet=0.0d0
+#endif
C
C Compute the side-chain and electrostatic interaction energy
C
else
esccor=0.0d0
endif
+#ifdef FOURBODY
C print *,"PRZED MULIt"
c print *,"Processor",myrank," computed Usccorr"
C
c & n_corr1
c call flush(iout)
endif
+#endif
c print *,"Processor",myrank," computed Ucorr"
c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode
if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
C print *,"przed lipidami"
if (wliptran.gt.0) then
call Eliptransfer(eliptran)
+ else
+ eliptran=0.0d0
endif
C print *,"za lipidami"
if (AFMlog.gt.0) then
end
c-------------------------------------------------------------------------------
subroutine sum_energy(energia,reduce)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifndef ISNAN
external proc_proc
#endif
#ifdef MPI
include "mpif.h"
+ integer ierr
+ double precision time00
#endif
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
include 'COMMON.TIME1'
logical reduce
+ integer i
+ double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
+ & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
+ & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6,
+ & eliptran,Eafmforce,Etube,
+ & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
+ double precision Uconst,etot
#ifdef MPI
if (nfgtasks.gt.1 .and. reduce) then
#ifdef DEBUG
end
c-------------------------------------------------------------------------------
subroutine sum_gradient
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifndef ISNAN
external proc_proc
#endif
#ifdef MPI
include 'mpif.h'
+ integer ierror,ierr
+ double precision time00,time01
#endif
double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres),
& glocbuf(4*maxres),gradbufc_sum(3,-1:maxres)
include 'COMMON.TIME1'
include 'COMMON.MAXGRAD'
include 'COMMON.SCCOR'
- include 'COMMON.MD'
+c include 'COMMON.MD'
+ include 'COMMON.QRESTR'
+ integer i,j,k
+ double precision scalar
+ double precision gvdwc_norm,gvdwc_scp_norm,gelc_norm,gvdwpp_norm,
+ &gradb_norm,ghpbc_norm,gradcorr_norm,gel_loc_norm,gcorr3_turn_norm,
+ &gcorr4_turn_norm,gradcorr5_norm,gradcorr6_norm,
+ &gcorr6_turn_norm,gsccorrc_norm,gscloc_norm,gvdwx_norm,
+ &gradx_scp_norm,ghpbx_norm,gradxorr_norm,gsccorrx_norm,
+ &gsclocx_norm
#ifdef TIMING
time01=MPI_Wtime()
#endif
gradcorr5_max=0.0d0
gradcorr6_max=0.0d0
gcorr6_turn_max=0.0d0
- gsccorc_max=0.0d0
+ gsccorrc_max=0.0d0
gscloc_max=0.0d0
gvdwx_max=0.0d0
gradx_scp_max=0.0d0
ghpbx_max=0.0d0
gradxorr_max=0.0d0
- gsccorx_max=0.0d0
+ gsccorrx_max=0.0d0
gsclocx_max=0.0d0
do i=1,nct
gvdwc_norm=dsqrt(scalar(gvdwc(1,i),gvdwc(1,i)))
if (gradcorr5_norm.gt.gradcorr5_max)
& gradcorr5_max=gradcorr5_norm
gradcorr6_norm=dsqrt(scalar(gradcorr6(1,i),gradcorr6(1,i)))
- if (gradcorr6_norm.gt.gradcorr6_max) gcorr6_max=gradcorr6_norm
+ if (gradcorr6_norm.gt.gradcorr6_max)gradcorr6_max=gradcorr6_norm
gcorr6_turn_norm=dsqrt(scalar(gcorr6_turn(1,i),
& gcorr6_turn(1,i)))
if (gcorr6_turn_norm.gt.gcorr6_turn_max)
& gcorr6_turn_max=gcorr6_turn_norm
- gsccorr_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i)))
- if (gsccorr_norm.gt.gsccorr_max) gsccorr_max=gsccorr_norm
+ gsccorrc_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i)))
+ if (gsccorrc_norm.gt.gsccorrc_max) gsccorrc_max=gsccorrc_norm
gscloc_norm=dsqrt(scalar(gscloc(1,i),gscloc(1,i)))
if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm
gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i)))
write (istat,'(1h#,21f10.2)') gvdwc_max,gvdwc_scp_max,
& gelc_max,gvdwpp_max,gradb_max,ghpbc_max,
& gradcorr_max,gel_loc_max,gcorr3_turn_max,gcorr4_turn_max,
- & gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorc_max,
+ & gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorrc_max,
& gscloc_max,gvdwx_max,gradx_scp_max,ghpbx_max,gradxorr_max,
- & gsccorx_max,gsclocx_max
+ & gsccorrx_max,gsclocx_max
close(istat)
if (gvdwc_max.gt.1.0d4) then
write (iout,*) "gvdwc gvdwx gradb gradbx"
end
c-------------------------------------------------------------------------------
subroutine rescale_weights(t_bath)
- implicit real*8 (a-h,o-z)
+ implicit none
+#ifdef MPI
+ include 'mpif.h'
+ integer ierror
+#endif
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.SBRIDGE'
include 'COMMON.CONTROL'
+ double precision t_bath
+ double precision facT,facT2,facT3,facT4,facT5
double precision kfac /2.4d0/
double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/
c facT=temp0/t_bath
end
C------------------------------------------------------------------------
subroutine enerprint(energia)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.SBRIDGE'
- include 'COMMON.MD'
+ include 'COMMON.QRESTR'
double precision energia(0:n_ene)
+ double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
+ & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
+ & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,
+ & eello_turn6,
+ & eliptran,Eafmforce,Etube,
+ & esaxs,ehomology_constr,edfator,edfanei,edfabet,etot
etot=energia(0)
evdw=energia(1)
evdw2=energia(2)
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
+#ifdef FOURBODY
& ecorr,wcorr,
- & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
- & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,
- & ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
+ & ecorr5,wcorr5,ecorr6,wcorr6,
+#endif
+ & eel_loc,wel_loc,eello_turn3,wturn3,
+ & eello_turn4,wturn4,
+#ifdef FOURBODY
+ & eello_turn6,wturn6,
+#endif
+ & esccor,wsccor,edihcnstr,
+ & ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforce,
& etube,wtube,esaxs,wsaxs,ehomology_constr,
& edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei,
& edfabet,wdfa_beta,
& 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
& 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6,
& ' (SS bridges & dist. cnstr.)'/
+#ifdef FOURBODY
& 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+#endif
& 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
& 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
& 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+#ifdef FOURBODY
& 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+#endif
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/
& 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
& estr,wbond,ebe,wang,
& escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
+#ifdef FOURBODY
& ecorr,wcorr,
- & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
- & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,
+ & ecorr5,wcorr5,ecorr6,wcorr6,
+#endif
+ & eel_loc,wel_loc,eello_turn3,wturn3,
+ & eello_turn4,wturn4,
+#ifdef FOURBODY
+ & eello_turn6,wturn6,
+#endif
+ & esccor,wsccor,edihcnstr,
& ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
& etube,wtube,esaxs,wsaxs,ehomology_constr,
& edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei,
& 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
& 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6,
& ' (SS bridges & dist. restr.)'/
+#ifdef FOURBODY
& 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+#endif
& 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
& 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
& 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+#ifdef FOURBODY
& 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+#endif
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/
& 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the LJ potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
+ double precision accur
include 'DIMENSIONS'
parameter (accur=1.0d-10)
include 'COMMON.GEO'
include 'COMMON.SBRIDGE'
include 'COMMON.NAMES'
include 'COMMON.IOUNITS'
+ include 'COMMON.SPLITELE'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
- dimension gg(3)
+ include 'COMMON.CONTMAT'
+#endif
+ double precision gg(3)
+ double precision evdw,evdwij
+ 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 fcont,fprimcont
+ double precision sscale,sscagrad
+ double precision boxshift
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do ikont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(ikont)
+ j=newcontlistj(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C Change 12/1/95
num_conti=0
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
+c do iint=1,nint_gr(i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
- do j=istart(i,iint),iend(i,iint)
+c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
C Change 12/1/95 to calculate four-body interactions
rij=xj*xj+yj*yj+zj*zj
rrij=1.0D0/rij
+ sqrij=dsqrt(rij)
+ sss1=sscale(sqrij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(sqrij,r_cut_int)
+
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj),
cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- evdw=evdw+evdwij
+ evdw=evdw+sss1*evdwij
C
C Calculate the components of the gradient in DC and X
C
- fac=-rrij*(e1+evdwij)
+ fac=-rrij*(e1+evdwij)*sss1
+ & +evdwij*sssgrad1/sqrij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
cgrad enddo
cgrad enddo
C
+#ifdef FOURBODY
C 12/1/95, revised on 5/20/97
C
C Calculate the contact function. The ith column of the array JCONT will
cd & i,j,(gacont(kk,num_conti,i),kk=1,3)
endif
endif
- enddo ! j
- enddo ! iint
+#endif
+c enddo ! j
+c enddo ! iint
C Change 12/1/95
+#ifdef FOURBODY
num_cont(i)=num_conti
+#endif
enddo ! i
do i=1,nct
do j=1,3
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the LJK potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
- dimension gg(3)
+ include 'COMMON.SPLITELE'
+ double precision gg(3)
+ double precision evdw,evdwij
+ integer i,j,k,itypi,itypj,itypi1,iint,ikont
+ 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
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do ikont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(ikont)
+ j=newcontlistj(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
r_inv_ij=dsqrt(rrij)
rij=1.0D0/r_inv_ij
+ sss1=sscale(rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(rij,r_cut_int)
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
C have you changed here?
cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- evdw=evdw+evdwij
+ evdw=evdw+evdwij*sss1
C
C Calculate the components of the gradient in DC and X
C
fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ & +evdwij*sssgrad1*r_inv_ij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
cgrad enddo
cgrad enddo
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
do i=1,nct
do j=1,3
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Berne-Pechukas potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SPLITELE'
+ integer icall
common /srutu/ icall
+ double precision evdw
+ 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 boxshift
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
lprn=.false.
c endif
ind=0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do ikont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(ikont)
+ j=newcontlistj(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
cd rrij=rrsave(ind)
cd endif
rij=dsqrt(rrij)
+ sss1=sscale(1.0d0/rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
C Calculate the angle-dependent terms of energy & contributions to derivatives.
call sc_angular
C Calculate whole angle-dependent part of epsilon and contributions
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij
+ evdw=evdw+sss1*evdwij
if (lprn) then
sigm=dabs(aa/bb)**(1.0D0/6.0D0)
epsi=bb**2/aa
fac=-expon*(e1+evdwij)
sigder=fac/sigsq
fac=rrij*fac
+ & +evdwij*sssgrad1/sss1*rij
C Calculate radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
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
- enddo ! j
- enddo ! iint
+! enddo ! j
+! enddo ! iint
enddo ! i
c stop
return
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Gay-Berne potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.SPLITELE'
include 'COMMON.SBRIDGE'
logical lprn
- integer xshift,yshift,zshift
-
+ double precision evdw
+ integer itypi,itypj,itypi1,iint,ind,ikont
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+ double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
+ & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
+ double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+ double precision boxshift
evdw=0.0D0
ccccc energy_dec=.false.
C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
C do xshift=-1,1
C do yshift=-1,1
C do zshift=-1,1
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do ikont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(ikont)
+ j=newcontlistj(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
-C Return atom into box, boxxsize is size of box in x dimension
-c 134 continue
-c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
-c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
-C Condition for being inside the proper box
-c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
-c & (xi.lt.((xshift-0.5d0)*boxxsize))) then
-c go to 134
-c endif
-c 135 continue
-c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
-c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
-C Condition for being inside the proper box
-c if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
-c & (yi.lt.((yshift-0.5d0)*boxysize))) then
-c go to 135
-c endif
-c 136 continue
-c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
-c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
-C Condition for being inside the proper box
-c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
-c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
-c go to 136
-c endif
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
C define scaling factor for lipids
C if (positi.le.0) positi=positi+boxzsize
C print *,i
C first for peptide groups
c for each residue check if it is in lipid or lipid water border area
- if ((zi.gt.bordlipbot)
- &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
- if (zi.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
-
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
C xi=xi+xshift*boxxsize
C yi=yi+yshift*boxysize
C zi=zi+zshift*boxzsize
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
c write(iout,*) "PRZED ZWYKLE", evdwij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
& 'evdw',i,j,evdwij,' ss'
C triple bond artifac removal
- do k=j+1,iend(i,iint)
+ do k=j+1,iend(i,iint)
C search over all next residues
- if (dyn_ss_mask(k)) then
+ if (dyn_ss_mask(k)) then
C check if they are cysteins
C write(iout,*) 'k=',k
c write(iout,*) "PRZED TRI", evdwij
- evdwij_przed_tri=evdwij
- call triple_ssbond_ene(i,j,k,evdwij)
+ evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
c if(evdwij_przed_tri.ne.evdwij) then
c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
c endif
c write(iout,*) "PO TRI", evdwij
C call the energy function that removes the artifical triple disulfide
C bond the soubroutine is located in ssMD.F
- evdw=evdw+evdwij
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
& 'evdw',i,j,evdwij,'tss'
- endif!dyn_ss_mask(k)
- enddo! k
+ endif!dyn_ss_mask(k)
+ enddo! k
ELSE
- ind=ind+1
- itypj=iabs(itype(j))
- if (itypj.eq.ntyp1) cycle
+ ind=ind+1
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
- dscj_inv=vbld_inv(j+nres)
+ dscj_inv=vbld_inv(j+nres)
c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
c & 1.0d0/vbld(j+nres)
c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
- sig0ij=sigma(itypi,itypj)
- chi1=chi(itypi,itypj)
- chi2=chi(itypj,itypi)
- chi12=chi1*chi2
- chip1=chip(itypi)
- chip2=chip(itypj)
- chip12=chip1*chip2
- alf1=alp(itypi)
- alf2=alp(itypj)
- alf12=0.5D0*(alf1+alf2)
+ sig0ij=sigma(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
C For diagnostics only!!!
c chi1=0.0D0
c chi2=0.0D0
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)
- yj=c(2,nres+j)
- zj=c(3,nres+j)
-C Return atom J into box the original box
-c 137 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 137
-c endif
-c 138 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 138
-c endif
-c 139 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 139
-c endif
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot)
- &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zj.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
- 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=c(1,nres+j)
+ 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
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 if (ssgradlipj.gt.0.0d0) print *,"??WTF??"
C print *,sslipi,sslipj,bordlipbot,zi,zj
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
- dxj=dc_norm(1,nres+j)
- dyj=dc_norm(2,nres+j)
- dzj=dc_norm(3,nres+j)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
C xj=xj-xi
C yj=yj-yi
C zj=zj-zi
c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
c write (iout,*) "j",j," dc_norm",
c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- rij=dsqrt(rrij)
- sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
- sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
-
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/rij,r_cut_int)
c write (iout,'(a7,4f8.3)')
c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
- if (sss.gt.0.0d0) then
+ if (sss.eq.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/rij,r_cut_int)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
- call sc_angular
- sigsq=1.0D0/sigsq
- sig=sig0ij*dsqrt(sigsq)
- rij_shift=1.0D0/rij-sig+sig0ij
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
c for diagnostics; uncomment
c rij_shift=1.2*sig0ij
C I hate to put IF's in the loops, but here don't have another choice!!!!
- if (rij_shift.le.0.0D0) then
- evdw=1.0D20
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
- return
- endif
- sigder=-sig*sigsq
+ return
+ endif
+ sigder=-sig*sigsq
c---------------------------------------------------------------
- rij_shift=1.0D0/rij_shift
- fac=rij_shift**expon
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
C here to start with
C if (c(i,3).gt.
- faclip=fac
- e1=fac*fac*aa
- e2=fac*bb
- evdwij=eps1*eps2rt*eps3rt*(e1+e2)
- eps2der=evdwij*eps3rt
- eps3der=evdwij*eps2rt
+ faclip=fac
+ e1=fac*fac*aa
+ e2=fac*bb
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij,
C &((sslipi+sslipj)/2.0d0+
C &(2.0d0-sslipi-sslipj)/2.0d0)
c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
- evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*sss
- if (lprn) then
- sigm=dabs(aa/bb)**(1.0D0/6.0D0)
- epsi=bb**2/aa
- write (iout,'(2(a3,i3,2x),17(0pf7.3))')
- & restyp(itypi),i,restyp(itypj),j,
- & epsi,sigm,chi1,chi2,chip1,chip2,
- & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
- & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
- & evdwij
- endif
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*sss
+ if (lprn) then
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+ & evdwij
+ endif
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'evdw',i,j,evdwij
+ if (energy_dec) write (iout,'(a,2i5,3f10.5)')
+ & 'r sss evdw',i,j,1.0d0/rij,sss,evdwij
C Calculate gradient components.
- e1=e1*eps1*eps2rt**2*eps3rt**2
- fac=-expon*(e1+evdwij)*rij_shift
- sigder=fac*sigder
- fac=rij*fac
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
c print '(2i4,6f8.4)',i,j,sss,sssgrad*
c & evdwij,fac,sigma(itypi,itypj),expon
- fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+ fac=fac+evdwij*sssgrad/sss*rij
c fac=0.0d0
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_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 gg_lipi(3)=0.0d0
C gg_lipj(3)=0.0d0
- gg(1)=xj*fac
- gg(2)=yj*fac
- gg(3)=zj*fac
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
C Calculate angular part of the gradient.
- call sc_grad
- endif
+c call sc_grad_scale(sss)
+ call sc_grad
ENDIF ! dyn_ss
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
C enddo ! zshift
C enddo ! yshift
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Gay-Berne-Vorobjev potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
- integer xshift,yshift,zshift
+ include 'COMMON.SPLITELE'
+ double precision boxshift
+ integer icall
common /srutu/ icall
logical lprn
+ double precision evdw
+ integer itypi,itypj,itypi1,iint,ind,ikont
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
+ & xi,yi,zi,fac_augm,e_augm
+ double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
+ & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip,sssgrad1
+ double precision dist,sscale,sscagrad,sscagradlip,sscalelip
evdw=0.0D0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
lprn=.false.
c if (icall.eq.0) lprn=.true.
ind=0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do ikont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(ikont)
+ j=newcontlistj(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
C define scaling factor for lipids
C if (positi.le.0) positi=positi+boxzsize
C print *,i
C first for peptide groups
c for each residue check if it is in lipid or lipid water border area
- if ((zi.gt.bordlipbot)
- &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
- if (zi.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
-
+ 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)
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
-C xj=c(1,nres+j)-xi
-C yj=c(2,nres+j)-yi
-C zj=c(3,nres+j)-zi
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot)
- &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zj.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
- 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=c(1,nres+j)
+ 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
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,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
- dxj=dc_norm(1,nres+j)
- dyj=dc_norm(2,nres+j)
- dzj=dc_norm(3,nres+j)
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- rij=dsqrt(rrij)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/rij,r_cut_int)
+ if (sss.eq.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/rij,r_cut_int)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
- call sc_angular
- sigsq=1.0D0/sigsq
- sig=sig0ij*dsqrt(sigsq)
- rij_shift=1.0D0/rij-sig+r0ij
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+r0ij
C I hate to put IF's in the loops, but here don't have another choice!!!!
- if (rij_shift.le.0.0D0) then
- evdw=1.0D20
- return
- endif
- sigder=-sig*sigsq
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
c---------------------------------------------------------------
- rij_shift=1.0D0/rij_shift
- fac=rij_shift**expon
- e1=fac*fac*aa
- e2=fac*bb
- evdwij=eps1*eps2rt*eps3rt*(e1+e2)
- eps2der=evdwij*eps3rt
- eps3der=evdwij*eps2rt
- fac_augm=rrij**expon
- e_augm=augm(itypi,itypj)*fac_augm
- evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij+e_augm
- if (lprn) then
- sigm=dabs(aa/bb)**(1.0D0/6.0D0)
- epsi=bb**2/aa
- write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa
+ e2=fac*bb
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij+e_augm
+ if (lprn) then
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
& chi1,chi2,chip1,chip2,
& eps1,eps2rt**2,eps3rt**2,
& om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
& evdwij+e_augm
- endif
+ endif
C Calculate gradient components.
- e1=e1*eps1*eps2rt**2*eps3rt**2
- fac=-expon*(e1+evdwij)*rij_shift
- sigder=fac*sigder
- fac=rij*fac-2*expon*rrij*e_augm
- fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ 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(1)=xj*fac
- gg(2)=yj*fac
- gg(3)=zj*fac
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
C Calculate angular part of the gradient.
- call sc_grad
- enddo ! j
- enddo ! iint
+c call sc_grad_scale(sss)
+ call sc_grad
+c enddo ! j
+c enddo ! iint
enddo ! i
end
C-----------------------------------------------------------------------------
include 'COMMON.SBRIDGE'
include 'COMMON.NAMES'
include 'COMMON.IOUNITS'
- include 'COMMON.CONTACTS'
+c include 'COMMON.CONTACTS'
dimension gg(3)
+ double precision boxshift
cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
evdw=0.0D0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do ikont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(ikont)
+ j=newcontlistj(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
+c do iint=1,nint_gr(i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
- do j=istart(i,iint),iend(i,iint)
+c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=boxshift(c(1,nres+j)-xi,boxxsize)
+ yj=boxshift(c(2,nres+j)-yi,boxysize)
+ zj=boxshift(c(3,nres+j)-zi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
r0ij=r0(itypi,itypj)
cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
cgrad enddo
cgrad enddo
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
return
end
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
- include 'COMMON.CONTACTS'
+c include 'COMMON.CONTACTS'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
dimension ggg(3)
- integer xshift,yshift,zshift
+ double precision boxshift
C write(iout,*) 'In EELEC_soft_sphere'
ees=0.0D0
evdw1=0.0D0
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
- sss=sscale(sqrt(rij))
- sssgrad=sscagrad(sqrt(rij))
+ sss=sscale(sqrt(rij),r_cut_int)
+ sssgrad=sscagrad(sqrt(rij),r_cut_int)
if (rij.lt.r0ijsq) then
evdw1ij=0.25d0*(rij-r0ijsq)**2
fac=rij-r0ijsq
endif
do j=1,2
do k=1,3
- do l=1,3
- uygrad(l,k,j,i)=vbld_inv_temp(j)*uygrad(l,k,j,i)
- uzgrad(l,k,j,i)=vbld_inv_temp(j)*uzgrad(l,k,j,i)
- enddo
- enddo
- enddo
- enddo
-#if defined(PARVEC) && defined(MPI)
- if (nfgtasks1.gt.1) then
- time00=MPI_Wtime()
-c print *,"Processor",fg_rank1,kolor1," ivec_start",ivec_start,
-c & " ivec_displ",(ivec_displ(i),i=0,nfgtasks1-1),
-c & " ivec_count",(ivec_count(i),i=0,nfgtasks1-1)
- call MPI_Allgatherv(uy(1,ivec_start),ivec_count(fg_rank1),
- & MPI_UYZ,uy(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ,
- & FG_COMM1,IERR)
- call MPI_Allgatherv(uz(1,ivec_start),ivec_count(fg_rank1),
- & MPI_UYZ,uz(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ,
- & FG_COMM1,IERR)
- call MPI_Allgatherv(uygrad(1,1,1,ivec_start),
- & ivec_count(fg_rank1),MPI_UYZGRAD,uygrad(1,1,1,1),ivec_count(0),
- & ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR)
- call MPI_Allgatherv(uzgrad(1,1,1,ivec_start),
- & ivec_count(fg_rank1),MPI_UYZGRAD,uzgrad(1,1,1,1),ivec_count(0),
- & ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR)
- time_gather=time_gather+MPI_Wtime()-time00
- endif
-#endif
-#ifdef DEBUG
- if (fg_rank.eq.0) then
- write (iout,*) "Arrays UY and UZ"
- do i=1,nres-1
- write (iout,'(i5,3f10.5,5x,3f10.5)') i,(uy(k,i),k=1,3),
- & (uz(k,i),k=1,3)
- enddo
- endif
-#endif
- return
- end
-C-----------------------------------------------------------------------------
- subroutine check_vecgrad
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.IOUNITS'
- include 'COMMON.GEO'
- include 'COMMON.VAR'
- include 'COMMON.LOCAL'
- include 'COMMON.CHAIN'
- include 'COMMON.VECTORS'
- dimension uygradt(3,3,2,maxres),uzgradt(3,3,2,maxres)
- dimension uyt(3,maxres),uzt(3,maxres)
- dimension uygradn(3,3,2),uzgradn(3,3,2),erij(3)
- double precision delta /1.0d-7/
- call vec_and_deriv
-cd do i=1,nres
-crc write(iout,'(2i5,2(3f10.5,5x))') i,1,dc_norm(:,i)
-crc write(iout,'(2i5,2(3f10.5,5x))') i,2,uy(:,i)
-crc write(iout,'(2i5,2(3f10.5,5x)/)')i,3,uz(:,i)
-cd write(iout,'(2i5,2(3f10.5,5x))') i,1,
-cd & (dc_norm(if90,i),if90=1,3)
-cd write(iout,'(2i5,2(3f10.5,5x))') i,2,(uy(if90,i),if90=1,3)
-cd write(iout,'(2i5,2(3f10.5,5x)/)')i,3,(uz(if90,i),if90=1,3)
-cd write(iout,'(a)')
-cd enddo
- do i=1,nres
- do j=1,2
- do k=1,3
- do l=1,3
- uygradt(l,k,j,i)=uygrad(l,k,j,i)
- uzgradt(l,k,j,i)=uzgrad(l,k,j,i)
- enddo
- enddo
- enddo
- enddo
- call vec_and_deriv
- do i=1,nres
- do j=1,3
- uyt(j,i)=uy(j,i)
- uzt(j,i)=uz(j,i)
- enddo
- enddo
- do i=1,nres
-cd write (iout,*) 'i=',i
- do k=1,3
- erij(k)=dc_norm(k,i)
- enddo
- do j=1,3
- do k=1,3
- dc_norm(k,i)=erij(k)
- enddo
- dc_norm(j,i)=dc_norm(j,i)+delta
-c fac=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
-c do k=1,3
-c dc_norm(k,i)=dc_norm(k,i)/fac
-c enddo
-c write (iout,*) (dc_norm(k,i),k=1,3)
-c write (iout,*) (erij(k),k=1,3)
- call vec_and_deriv
- do k=1,3
- uygradn(k,j,1)=(uy(k,i)-uyt(k,i))/delta
- uygradn(k,j,2)=(uy(k,i-1)-uyt(k,i-1))/delta
- uzgradn(k,j,1)=(uz(k,i)-uzt(k,i))/delta
- uzgradn(k,j,2)=(uz(k,i-1)-uzt(k,i-1))/delta
- enddo
-c write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)')
-c & j,(uzgradt(k,j,1,i),k=1,3),(uzgradn(k,j,1),k=1,3),
-c & (uzgradt(k,j,2,i-1),k=1,3),(uzgradn(k,j,2),k=1,3)
- enddo
- do k=1,3
- dc_norm(k,i)=erij(k)
+ do l=1,3
+ uygrad(l,k,j,i)=vbld_inv_temp(j)*uygrad(l,k,j,i)
+ uzgrad(l,k,j,i)=vbld_inv_temp(j)*uzgrad(l,k,j,i)
+ enddo
+ enddo
enddo
-cd do k=1,3
-cd write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)')
-cd & k,(uygradt(k,l,1,i),l=1,3),(uygradn(k,l,1),l=1,3),
-cd & (uygradt(k,l,2,i-1),l=1,3),(uygradn(k,l,2),l=1,3)
-cd write (iout,'(i5,3f8.5,3x,3f8.5,5x,3f8.5,3x,3f8.5)')
-cd & k,(uzgradt(k,l,1,i),l=1,3),(uzgradn(k,l,1),l=1,3),
-cd & (uzgradt(k,l,2,i-1),l=1,3),(uzgradn(k,l,2),l=1,3)
-cd write (iout,'(a)')
-cd enddo
enddo
+#if defined(PARVEC) && defined(MPI)
+ if (nfgtasks1.gt.1) then
+ time00=MPI_Wtime()
+c print *,"Processor",fg_rank1,kolor1," ivec_start",ivec_start,
+c & " ivec_displ",(ivec_displ(i),i=0,nfgtasks1-1),
+c & " ivec_count",(ivec_count(i),i=0,nfgtasks1-1)
+ call MPI_Allgatherv(uy(1,ivec_start),ivec_count(fg_rank1),
+ & MPI_UYZ,uy(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ,
+ & FG_COMM1,IERR)
+ call MPI_Allgatherv(uz(1,ivec_start),ivec_count(fg_rank1),
+ & MPI_UYZ,uz(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ,
+ & FG_COMM1,IERR)
+ call MPI_Allgatherv(uygrad(1,1,1,ivec_start),
+ & ivec_count(fg_rank1),MPI_UYZGRAD,uygrad(1,1,1,1),ivec_count(0),
+ & ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR)
+ call MPI_Allgatherv(uzgrad(1,1,1,ivec_start),
+ & ivec_count(fg_rank1),MPI_UYZGRAD,uzgrad(1,1,1,1),ivec_count(0),
+ & ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR)
+ time_gather=time_gather+MPI_Wtime()-time00
+ endif
+#endif
+#ifdef DEBUG
+ if (fg_rank.eq.0) then
+ write (iout,*) "Arrays UY and UZ"
+ do i=1,nres-1
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(uy(k,i),k=1,3),
+ & (uz(k,i),k=1,3)
+ enddo
+ endif
+#endif
return
end
C--------------------------------------------------------------------------
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
- include 'COMMON.CONTACTS'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
#else
do i=3,nres+1
#endif
- if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ ii=ireschain(i-2)
+c write (iout,*) "i",i,i-2," ii",ii
+ if (ii.eq.0) cycle
+ innt=chain_border(1,ii)
+ inct=chain_border(2,ii)
+c write (iout,*) "i",i,i-2," ii",ii," innt",innt," inct",inct
+c if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (i.gt. innt+2 .and. i.lt.inct+2) then
iti = itype2loc(itype(i-2))
else
iti=nloctyp
endif
c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
- if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ if (i.gt. innt+1 .and. i.lt.inct+1) then
iti1 = itype2loc(itype(i-1))
else
iti1=nloctyp
endif
-c write(iout,*),i
+c write(iout,*),"i",i,i-2," iti",itype(i-2),iti,
+c & " iti1",itype(i-1),iti1
#ifdef NEWCORR
cost1=dcos(theta(i-1))
sint1=dsin(theta(i-1))
write (iout,*) 'theta=', theta(i-1)
#endif
#else
- if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (i.gt. innt+2 .and. i.lt.inct+2) then
+c if (i.gt. nnt+2 .and. i.lt.nct+2) then
iti = itype2loc(itype(i-2))
else
iti=nloctyp
write(iout,*) 'b2=',(b2(k,i-2),k=1,2)
#endif
enddo
+ mu=0.0d0
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do i=3,nres+1
#endif
- if (i .lt. nres+1) then
+c if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
+ if (i .lt. nres+1 .and. itype(i-1).lt.ntyp1) then
sin1=dsin(phi(i))
cos1=dcos(phi(i))
sintab(i-2)=sin1
Ug2(2,1,i-2)=0.0d0
Ug2(2,2,i-2)=0.0d0
endif
- if (i .gt. 3 .and. i .lt. nres+1) then
+ if (i .gt. 3) then
obrot_der(1,i-2)=-sin1
obrot_der(2,i-2)= cos1
Ugder(1,1,i-2)= sin1
Ug2der(2,2,i-2)=0.0d0
endif
c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
- if (i.gt. nnt+2 .and. i.lt.nct+2) then
+c if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (i.gt.nnt+2 .and.i.lt.nct+2) then
iti = itype2loc(itype(i-2))
else
iti=nloctyp
c write(iout,*) "Macierz EUG",
c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2),
c & eug(2,2,i-2)
+#ifdef FOURBODY
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
& then
call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2))
call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2))
call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2))
endif
+#endif
else
do k=1,2
Ub2(k,i-2)=0.0d0
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
cd write (iout,*) 'mu',i-2,mu(:,i-2)
+#ifdef FOURBODY
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
& then
call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
call matmat2(EUg(1,1,i-2),DD(1,1,i-1),EUgD(1,1,i-2))
call matmat2(EUgder(1,1,i-2),DD(1,1,i-1),EUgDder(1,1,i-2))
endif
+#endif
enddo
+#ifdef FOURBODY
C Matrices dependent on two consecutive virtual-bond dihedrals.
C The order of matrices is from left to right.
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
call matmat2(auxmat(1,1),EUg(1,1,i),Ug2DtEUgder(1,1,1,i))
enddo
endif
+#endif
#if defined(MPI) && defined(PARMAT)
#ifdef DEBUG
c if (fg_rank.eq.0) then
call MPI_Allgatherv(sintab2(ivec_start),ivec_count(fg_rank1),
& MPI_DOUBLE_PRECISION,sintab2(1),ivec_count(0),ivec_displ(0),
& MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
+#ifdef FOURBODY
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
& then
call MPI_Allgatherv(Ctobr(1,ivec_start),ivec_count(fg_rank1),
& MPI_MAT2,Ug2DtEUgder(1,1,1,1),ivec_count(0),ivec_displ(0),
& MPI_MAT2,FG_COMM1,IERR)
endif
+#endif
#else
c Passes matrix info through the ring
isend=fg_rank1
& iprev,6600+irecv,FG_COMM,status,IERR)
c write (iout,*) "Gather PRECOMP12"
c call flush(iout)
+#ifdef FOURBODY
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
& then
call MPI_SENDRECV(ug2db1t(1,ivec_displ(isend)+1),1,
& Ug2DtEUgder(1,1,1,ivec_displ(irecv)+1),1,
& MPI_PRECOMP23(lenrecv),
& iprev,9900+irecv,FG_COMM,status,IERR)
+#endif
c write (iout,*) "Gather PRECOMP23"
c call flush(iout)
endif
cd enddo
return
end
-C--------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
C
C This subroutine calculates the average interaction energy and its gradient
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+#endif
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
eello_turn3=0.0d0
eello_turn4=0.0d0
ind=0
+#ifdef FOURBODY
do i=1,nres
num_cont_hb(i)=0
enddo
+#endif
cd print '(a)','Enter EELEC'
cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
do i=1,nres
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
num_conti=0
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo
do i=iturn4_start,iturn4_end
if (i.lt.1) cycle
c & (zmedi.lt.((-0.5d0)*boxzsize))) then
c go to 196
c endif
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
-
+ call to_box(xmedi,ymedi,zmedi)
+#ifdef FOURBODY
num_conti=num_cont_hb(i)
+#endif
c write(iout,*) "JESTEM W PETLI"
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
& call eturn4(i,eello_turn4)
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo ! i
C Loop over all neighbouring boxes
C do xshift=-1,1
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
CTU KURWA
- do i=iatel_s,iatel_e
+c do i=iatel_s,iatel_e
+ do ikont=g_listpp_start,g_listpp_end
+ i=newcontlistppi(ikont)
+ j=newcontlistppj(ikont)
C do i=75,75
c if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
C xmedi=xmedi+xshift*boxxsize
C ymedi=ymedi+yshift*boxysize
C zmedi=zmedi+zshift*boxzsize
c endif
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+#ifdef FOURBODY
num_conti=num_cont_hb(i)
+#endif
C I TU KURWA
- do j=ielstart(i),ielend(i)
+c do j=ielstart(i),ielend(i)
C do j=16,17
C write (iout,*) i,j
C if (j.le.1) cycle
- if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+ if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
c & .or.((j+2).gt.nres)
c & .or.((j-1).le.0)
c & .or.itype(j+2).eq.ntyp1
c & .or.itype(j-1).eq.ntyp1
&) cycle
- call eelecij(i,j,ees,evdw1,eel_loc)
- enddo ! j
+ call eelecij(i,j,ees,evdw1,eel_loc)
+c enddo ! j
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo ! i
C enddo ! zshift
C enddo ! yshift
end
C-------------------------------------------------------------------------------
subroutine eelecij(i,j,ees,evdw1,eel_loc)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include "mpif.h"
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+#endif
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
include 'COMMON.SPLITELE'
include 'COMMON.SHIELD'
- dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+ double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
& aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
& gmuij2(4),gmuji2(4)
+ double precision dxi,dyi,dzi
+ double precision dx_normi,dy_normi,dz_normi,aux
+ integer j1,j2,lll,num_conti
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
+ integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ilist,iresshield
+ double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
+ double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
+ double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
+ & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
+ & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
+ & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
+ & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
+ & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
+ & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
+ & ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield
+ double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji
+ double precision xmedi,ymedi,zmedi
+ double precision sscale,sscagrad,scalar
+ double precision boxshift
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
double precision scal_el /1.0d0/
double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
& 0.0d0,1.0d0,0.0d0,
& 0.0d0,0.0d0,1.0d0/
- integer xshift,yshift,zshift
c time00=MPI_Wtime()
cd write (iout,*) "eelecij",i,j
c ind=ind+1
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
c 174 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 174
-c endif
-c 175 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 175
-c endif
-c 176 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 176
-c endif
-C endif !endPBC condintion
-C xj=xj-xmedi
-C yj=yj-ymedi
-C zj=zj-zmedi
rij=xj*xj+yj*yj+zj*zj
- sss=sscale(sqrt(rij))
- sssgrad=sscagrad(sqrt(rij))
+ sss=sscale(dsqrt(rij),r_cut_int)
+ if (sss.eq.0.0d0) return
+ sssgrad=sscagrad(dsqrt(rij),r_cut_int)
c if (sss.gt.0.0d0) then
rrmij=1.0D0/rij
rij=dsqrt(rij)
fac_shield(i)=1.0
fac_shield(j)=1.0
eesij=(el1+el2)
- ees=ees+eesij
+ ees=ees+eesij*sss
endif
evdw1=evdw1+evdwij*sss
cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
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,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)
endif
C
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
- ggg(1)=facel*xj
- ggg(2)=facel*yj
- ggg(3)=facel*zj
+ aux=facel*sss+rmij*sssgrad*eesij
+ 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
iresshield=shield_list(ilist,j)
do k=1,3
rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
- & *2.0
+ & *2.0*sss
gshieldx(k,iresshield)=gshieldx(k,iresshield)+
& rlocshield
- & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss
gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
do k=1,3
gshieldc(k,i)=gshieldc(k,i)+
- & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
gshieldc(k,j)=gshieldc(k,j)+
- & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
gshieldc(k,i-1)=gshieldc(k,i-1)+
- & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
gshieldc(k,j-1)=gshieldc(k,j-1)+
- & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
enddo
endif
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
- if (sss.gt.0.0) then
- ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
- ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
- ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
- else
- ggg(1)=0.0
- ggg(2)=0.0
- ggg(3)=0.0
- endif
+ facvdw=facvdw+sssgrad*rmij*evdwij
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
cgrad enddo
#else
C MARYSIA
- facvdw=(ev1+evdwij)*sss
+ facvdw=(ev1+evdwij)
facel=(el1+eesij)
fac1=fac
- fac=-3*rrmij*(facvdw+facvdw+facel)
+ fac=-3*rrmij*(facvdw+facvdw+facel)*sss
+ & +(evdwij+eesij)*sssgrad*rrmij
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
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
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
do k=1,3
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))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss
& *fac_shield(i)**2*fac_shield(j)**2
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))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss
& *fac_shield(i)**2*fac_shield(j)**2
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)
+ & *fac_shield(i)*fac_shield(j)*sss
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)
+ & *fac_shield(i)*fac_shield(j)*sss
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)
+ & *fac_shield(i)*fac_shield(j)*sss
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
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)
+ & *fac_shield(i)*fac_shield(j)*sss
#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
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
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
+ ggg(2)=aux*yj
+ ggg(3)=aux*zj
do l=1,3
- ggg(l)=(agg(l,1)*muij(1)+
+ 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)
+ & *fac_shield(i)*fac_shield(j)*sss
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)
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
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
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
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
enddo
ENDIF
C Change 12/26/95 to calculate four-body contributions to H-bonding energy
c if (j.gt.i+1 .and. num_conti.le.maxconts) then
+#ifdef FOURBODY
if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
& .and. num_conti.le.maxconts) then
c write (iout,*) i,j," entered corr"
C fac_shield(j)=0.6d0
endif
ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
C Diagnostics. Comment out or remove after debugging!
c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
enddo
gggp(1)=gggp(1)+ees0pijp*xj
+ & +ees0p(num_conti,i)/sss*rmij*xj*sssgrad
gggp(2)=gggp(2)+ees0pijp*yj
+ & +ees0p(num_conti,i)/sss*rmij*yj*sssgrad
gggp(3)=gggp(3)+ees0pijp*zj
+ & +ees0p(num_conti,i)/sss*rmij*zj*sssgrad
gggm(1)=gggm(1)+ees0mijp*xj
+ & +ees0m(num_conti,i)/sss*rmij*xj*sssgrad
gggm(2)=gggm(2)+ees0mijp*yj
+ & +ees0m(num_conti,i)/sss*rmij*yj*sssgrad
gggm(3)=gggm(3)+ees0mijp*zj
+ & +ees0m(num_conti,i)/sss*rmij*zj*sssgrad
C Derivatives due to the contact function
gacont_hbr(1,num_conti,i)=fprimcont*xj
gacont_hbr(2,num_conti,i)=fprimcont*yj
gacontp_hb1(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontp_hb2(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontp_hb3(k,num_conti,i)=gggp(k)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb1(k,num_conti,i)=!ghalfm
& +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb2(k,num_conti,i)=!ghalfm
& +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb3(k,num_conti,i)=gggm(k)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
enddo
C Diagnostics. Comment out or remove after debugging!
endif ! num_conti.le.maxconts
endif ! fcont.gt.0
endif ! j.gt.i+1
+#endif
if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
do k=1,4
do l=1,3
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
- include 'COMMON.CONTACTS'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
- include 'COMMON.CONTACTS'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
dimension ggg(3)
- integer xshift,yshift,zshift
+ double precision boxshift
evdw2=0.0D0
evdw2_14=0.0d0
r0_scp=4.5d0
C do xshift=-1,1
C do yshift=-1,1
C do zshift=-1,1
- do i=iatscp_s,iatscp_e
+c do i=iatscp_s,iatscp_e
+ do ikont=g_listscp_start,g_listscp_end
+ i=newcontlistscpi(ikont)
+ j=newcontlistscpj(ikont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
c go to 136
c endif
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
C xi=xi+xshift*boxxsize
C yi=yi+yshift*boxysize
C zi=zi+zshift*boxzsize
- do iint=1,nscp_gr(i)
+c do iint=1,nscp_gr(i)
- do j=iscpstart(i,iint),iscpend(i,iint)
+c do j=iscpstart(i,iint),iscpend(i,iint)
if (itype(j).eq.ntyp1) cycle
itypj=iabs(itype(j))
C Uncomment following three lines for SC-p interactions
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
-c 174 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 174
-c endif
-c 175 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-cC Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 175
-c endif
-c 176 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 176
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
-c c endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
C xj=xj-xi
C yj=yj-yi
C zj=zj-zi
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
-cgrad if (j.lt.i) then
-cd write (iout,*) 'j<i'
-C Uncomment following three lines for SC-p interactions
-c do k=1,3
-c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
-c enddo
-cgrad else
-cd write (iout,*) 'j>i'
-cgrad do k=1,3
-cgrad ggg(k)=-ggg(k)
-C Uncomment following line for SC-p interactions
-c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
-cgrad enddo
-cgrad endif
-cgrad do k=1,3
-cgrad gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
-cgrad enddo
-cgrad kstart=min0(i+1,j)
-cgrad kend=max0(i-1,j-1)
-cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
-cd write (iout,*) ggg(1),ggg(2),ggg(3)
-cgrad do k=kstart,kend
-cgrad do l=1,3
-cgrad gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
-cgrad enddo
-cgrad enddo
do k=1,3
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- enddo
+c enddo
- enddo ! iint
+c enddo ! iint
enddo ! i
C enddo !zshift
C enddo !yshift
C peptide-group centers and side chains and its gradient in virtual-bond and
C side-chain vectors.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
- integer xshift,yshift,zshift
- dimension ggg(3)
+ double precision ggg(3)
+ integer i,iint,j,k,iteli,itypj,subchap,ikont
+ double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
+ & fac,e1,e2,rij
+ double precision evdw2,evdw2_14,evdwij
+ double precision sscale,sscagrad
+ double precision boxshift
evdw2=0.0D0
evdw2_14=0.0d0
c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
C do xshift=-1,1
C do yshift=-1,1
C do zshift=-1,1
- if (energy_dec) write (iout,*) "escp:",r_cut,rlamb
- do i=iatscp_s,iatscp_e
+ if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
+c do i=iatscp_s,iatscp_e
+ do ikont=g_listscp_start,g_listscp_end
+ i=newcontlistscpi(ikont)
+ j=newcontlistscpj(ikont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
-c xi=xi+xshift*boxxsize
-c yi=yi+yshift*boxysize
-c zi=zi+zshift*boxzsize
-c print *,xi,yi,zi,'polozenie i'
-C Return atom into box, boxxsize is size of box in x dimension
-c 134 continue
-c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
-c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
-C Condition for being inside the proper box
-c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
-c & (xi.lt.((xshift-0.5d0)*boxxsize))) then
-c go to 134
-c endif
-c 135 continue
-c print *,xi,boxxsize,"pierwszy"
-
-c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
-c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
-C Condition for being inside the proper box
-c if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
-c & (yi.lt.((yshift-0.5d0)*boxysize))) then
-c go to 135
-c endif
-c 136 continue
-c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
-c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
-C Condition for being inside the proper box
-c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
-c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
-c go to 136
-c endif
- do iint=1,nscp_gr(i)
+ call to_box(xi,yi,zi)
+c do iint=1,nscp_gr(i)
- do j=iscpstart(i,iint),iscpend(i,iint)
+c do j=iscpstart(i,iint),iscpend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
C Uncomment following three lines for SC-p interactions
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
-c 174 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 174
-c endif
-c 175 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-cC Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 175
-c endif
-c 176 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 176
-c endif
-CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
c print *,xj,yj,zj,'polozenie j'
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
c print *,rrij
- sss=sscale(1.0d0/(dsqrt(rrij)))
+ sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz'
c if (sss.eq.0) print *,'czasem jest OK'
if (sss.le.0.0d0) cycle
- sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
+ sssgrad=sscagrad(1.0d0/(dsqrt(rrij)),r_cut_int)
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
endif
evdwij=e1+e2
evdw2=evdw2+evdwij*sss
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
- & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
+ if (energy_dec) write (iout,'(a6,2i5,3f7.3,2i3,3e11.3)')
+ & 'evdw2',i,j,1.0d0/dsqrt(rrij),sss,
+ & evdwij,iteli,itypj,fac,aad(itypj,iteli),
& bad(itypj,iteli)
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
c endif !endif for sscale cutoff
- enddo ! j
+c enddo ! j
- enddo ! iint
+c enddo ! iint
enddo ! i
c enddo !zshift
c enddo !yshift
estr=0.0d0
estr1=0.0d0
do i=ibondp_start,ibondp_end
+c 3/4/2020 Adam: removed dummy bond graient if Calpha and SC coords are
+c used
+#ifdef FIVEDIAG
+ if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle
+ diff = vbld(i)-vbldp0
+#else
if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
c estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
c do j=1,3
c & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
c else
C Checking if it involves dummy (NH3+ or COO-) group
- if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+ if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
C YES vbldpDUM is the equlibrium length of spring for Dummy atom
- diff = vbld(i)-vbldpDUM
- if (energy_dec) write(iout,*) "dum_bond",i,diff
- else
-C NO vbldp0 is the equlibrium lenght of spring for peptide group
- diff = vbld(i)-vbldp0
- endif
- if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+ diff = vbld(i)-vbldpDUM
+ if (energy_dec) write(iout,*) "dum_bond",i,diff
+ else
+C NO vbldp0 is the equlibrium length of spring for peptide group
+ diff = vbld(i)-vbldp0
+ endif
+#endif
+ if (energy_dec) write (iout,'(a7,i5,4f7.3)')
& "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
do j=1,3
& phii1*rad2deg,ethetai
c lprn1=.false.
etheta=etheta+ethetai
- if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
- & 'ebend',i,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
c & dscp1,dscp2,sumene
c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
- if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
- & 'escloc',i,sumene
-c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
+ if (energy_dec) write (2,*) "i",i," itype",itype(i)," it",it,
+ & " escloc",sumene,escloc,it,itype(i)
c & ,zz,xx,yy
c#define DEBUG
#ifdef DEBUG
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.TORCNSTR'
- include 'COMMON.CONTROL'
logical lprn
C Set lprn=.true. for debugging
lprn=.false.
& (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
& (itype(i+1).eq.ntyp1)) cycle
C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
- etors_d_ii=0.0D0
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
itori2=itortyp(itype(i))
sinphi2=dsin(j*phii1)
etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
& v2cij*cosphi2+v2sij*sinphi2
- if (energy_dec) etors_d_ii=etors_d_ii+
- & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
sinphi1m2=dsin(l*phii-(k-l)*phii1)
etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
& v1sdij*sinphi1p2+v2sdij*sinphi1m2
- if (energy_dec) etors_d_ii=etors_d_ii+
- & v1cdij*cosphi1p2+v2cdij*cosphi1m2+
- & v1sdij*sinphi1p2+v2sdij*sinphi1m2
gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
& -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
& -v1cdij*sinphi1p2+v2cdij*sinphi1m2)
enddo
enddo
- if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
- & 'etor_d',i,etors_d_ii
gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
enddo
c----------------------------------------------------------------------------
c MODELLER restraint function
subroutine e_modeller(ehomology_constr)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
- integer nnn, i, j, k, ki, irec, l
+ double precision ehomology_constr
+ integer nnn,i,ii,j,k,ijk,jik,ki,kk,nexl,irec,l
integer katy, odleglosci, test7
real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
real*8 Eval,Erot
double precision, dimension (max_template) ::
& gtheta,dscdiff,uscdiffk,guscdiff2,guscdiff3,
& theta_diff
+ double precision sum_godl,sgodl,grad_odl3,ggodl,sum_gdih,
+ & sum_guscdiff,sum_sgdih,sgdih,grad_dih3,usc_diff_i,dxx,dyy,dzz,
+ & betai,sum_sgodl,dij
+ double precision dist,pinorm
c
-
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.GEO'
include 'COMMON.INTERACT'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
- include 'COMMON.MD'
+c include 'COMMON.MD'
include 'COMMON.CONTROL'
+ include 'COMMON.HOMOLOGY'
+ include 'COMMON.QRESTR'
c
c From subroutine Econstr_back
c
esccor=0.0D0
do i=itau_start,itau_end
if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
+ esccor_ii=0.0D0
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
c write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
phii=phi(i)
do intertyp=1,3 !intertyp
- esccor_ii=0.0D0
cc Added 09 May 2012 (Adasko)
cc Intertyp means interaction type of backbone mainchain correlation:
c 1 = SC...Ca...Ca...Ca
v2ij=v2sccor(j,intertyp,isccori,isccori1)
cosphi=dcos(j*tauangle(intertyp,i))
sinphi=dsin(j*tauangle(intertyp,i))
- if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
- if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)')
- & 'esccor',i,intertyp,esccor_ii
c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
if (lprn)
return
end
+#ifdef FOURBODY
c----------------------------------------------------------------------------
subroutine multibody(ecorr)
C This subroutine calculates multi-body contributions to energy following
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
double precision gx(3),gx1(3)
logical lprn
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.SHIELD'
double precision gx(3),gx1(3)
logical lprn
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.CONTROL'
include 'COMMON.LOCAL'
double precision gx(3),gx1(3),time00
parameter (max_cont=maxconts)
parameter (max_dim=26)
include "COMMON.CONTACTS"
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
double precision zapas(max_dim,maxconts,max_fg_procs),
& zapas_recv(max_dim,maxconts,max_fg_procs)
common /przechowalnia/ zapas
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.CHAIN'
include 'COMMON.CONTROL'
include 'COMMON.SHIELD'
parameter (max_cont=maxconts)
parameter (max_dim=70)
include "COMMON.CONTACTS"
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
double precision zapas(max_dim,maxconts,max_fg_procs),
& zapas_recv(max_dim,maxconts,max_fg_procs)
common /przechowalnia/ zapas
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.SHIELD'
include 'COMMON.CONTROL'
double precision gx(3),gx1(3)
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
cd write (2,*) 'eel_turn6',ekont*eel_turn6
return
end
-
C-----------------------------------------------------------------------------
+#endif
double precision function scalar(u,v)
!DIR$ INLINEALWAYS scalar
#ifndef OSF
include 'COMMON.INTERACT'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
- include 'COMMON.MD'
+c include 'COMMON.MD'
+#ifdef LANG0
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+#else
+ include 'COMMON.LANGEVIN'
+#endif
include 'COMMON.CONTROL'
+ include 'COMMON.SAXS'
include 'COMMON.NAMES'
include 'COMMON.TIME1'
include 'COMMON.FFIELD'
include 'COMMON.INTERACT'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
- include 'COMMON.MD'
+c include 'COMMON.MD'
+#ifdef LANG0
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+#else
+ include 'COMMON.LANGEVIN'
+#endif
include 'COMMON.CONTROL'
+ include 'COMMON.SAXS'
include 'COMMON.NAMES'
include 'COMMON.TIME1'
include 'COMMON.FFIELD'
endif
return
end
+c------------------------------------------------------------------------
+ double precision function boxshift(x,boxsize)
+ implicit none
+ double precision x,boxsize
+ double precision xtemp
+ xtemp=dmod(x,boxsize)
+ if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+ boxshift=xtemp-boxsize
+ else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+ boxshift=xtemp+boxsize
+ else
+ boxshift=xtemp
+ endif
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine closest_img(xi,yi,zi,xj,yj,zj)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ integer xshift,yshift,zshift,subchap
+ double precision dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine to_box(xi,yi,zi)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ double precision xi,yi,zi
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0.0d0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0.0d0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0.0d0) zi=zi+boxzsize
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ double precision xi,yi,zi,sslipi,ssgradlipi
+ double precision fracinbuf
+ double precision sscalelip,sscagradlip
+
+ if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+ return
+ end
projects / unres.git / blobdiff