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
+ double precision time01
#ifdef MPI
c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
c & " nfgtasks",nfgtasks
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
c call chainbuild_cart
endif
-#ifndef DFA
- edfadis=0.0d0
- edfator=0.0d0
- edfanei=0.0d0
- edfabet=0.0d0
+ 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
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
+ call make_SCp_inter_list
+c write (iout,*) "Finished make_SCp_inter_list"
+c call flush(iout)
+ call make_SCSC_inter_list
+c write (iout,*) "Finished make_SCSC_inter_list"
+c call flush(iout)
+ call make_pp_inter_list
+c write (iout,*) "Finished make_pp_inter_list"
+c call flush(iout)
+c call make_pp_vdw_inter_list
+c write (iout,*) "Finished make_pp_vdw_inter_list"
+c call flush(iout)
+#ifdef TIMING_ENE
+ time_list=time_list+MPI_Wtime()-time01
#endif
+ 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
C print *,ipot
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
goto (101,102,103,104,105,106) ipot
C Lennard-Jones potential.
101 call elj(evdw)
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+#ifdef TIMING_ENE
+ time_evdw=time_evdw+MPI_Wtime()-time01
+#endif
#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
time_vec=time_vec+MPI_Wtime()-time01
#endif
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
C Introduction of shielding effect first for each peptide group
C the shielding factor is set this factor is describing how each
C peptide group is shielded by side-chains
c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
c & eello_turn4)
endif
+#ifdef TIMING_ENE
+ time_eelec=time_eelec+MPI_Wtime()-time01
+#endif
c#ifdef TIMING
c time_enecalc=time_enecalc+MPI_Wtime()-time00
c#endif
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
C
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
if (ipot.lt.6) then
if(wscp.gt.0d0) then
call escp(evdw2,evdw2_14)
c write (iout,*) "Soft-sphere SCP potential"
call escp_soft_sphere(evdw2,evdw2_14)
endif
+#ifdef TIMING_ENE
+ time_escp=time_escp+MPI_Wtime()-time01
+#endif
c
c Calculate the bond-stretching energy
c
c write (iout,*) "MULTIBODY_HB ecorr",ecorr,ecorr5,ecorr6,n_corr,
c & n_corr1
c call flush(iout)
+ else
+ ecorr=0.0d0
+ ecorr5=0.0d0
+ ecorr6=0.0d0
+ eturn6=0.0d0
endif
+#else
+ ecorr=0.0d0
+ ecorr5=0.0d0
+ ecorr6=0.0d0
+ eturn6=0.0d0
#endif
c print *,"Processor",myrank," computed Ucorr"
c write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode
call AFMforce(Eafmforce)
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
+ else
+ Eafmforce=0.0d0
endif
if (TUBElog.eq.1) then
C print *,"just before call"
call calctube(Etube)
- elseif (TUBElog.eq.2) then
+ elseif (TUBElog.eq.2) then
call calctube2(Etube)
- else
- Etube=0.0d0
- endif
+ else
+ Etube=0.0d0
+ endif
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#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
enddo
-#ifdef DEBUG
- write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
- write (iout,*) (i," jgrad_start",jgrad_start(i),
- & " jgrad_end ",jgrad_end(i),
- & i=igrad_start,igrad_end)
-#endif
+c#ifdef DEBUG
+c write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
+c write (iout,*) (i," jgrad_start",jgrad_start(i),
+c & " jgrad_end ",jgrad_end(i),
+c & i=igrad_start,igrad_end)
+c#endif
c
c Obsolete and inefficient code; we can make the effort O(n) and, therefore,
c do not parallelize this part.
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
#endif
double precision gg(3)
double precision evdw,evdwij
- integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+ 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 fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
+ & faclip
+ double precision sscale,sscagrad,sscagradlip,sscalelip
+ double precision gg_lipi(3),gg_lipj(3)
+ double precision boxshift
+ external boxshift
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
- do i=iatsc_s,iatsc_e
+ 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)
+ 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)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
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)
+ 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 Change 12/1/95 to calculate four-body interactions
rij=xj*xj+yj*yj+zj*zj
rrij=1.0D0/rij
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
endif
endif
#endif
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
C Change 12/1/95
#ifdef FOURBODY
num_cont(i)=num_conti
include 'COMMON.SPLITELE'
double precision gg(3)
double precision evdw,evdwij
- integer i,j,k,itypi,itypj,itypi1,iint
+ 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
+ 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
- do i=iatsc_s,iatsc_e
+ 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)
+ 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)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
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)
+ 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/(xj*xj+yj*yj+zj*zj)
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
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
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
integer icall
common /srutu/ icall
double precision evdw
- integer itypi,itypj,itypi1,iint,ind
+ 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
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)
+ 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
- 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)
+ 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)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
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
- enddo ! j
- enddo ! iint
+! enddo ! j
+! enddo ! iint
enddo ! i
c stop
return
include 'COMMON.SPLITELE'
include 'COMMON.SBRIDGE'
logical lprn
- integer xshift,yshift,zshift,subchap
double precision evdw
- integer itypi,itypj,itypi1,iint,ind
+ 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,dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+ double precision boxshift
+ double precision x0,y0,r012,rij12,facx0,
+ & facx02,afacx0,bfacx0,abfacx0,Afac,BBfac,Afacsig,Bfacsig
+c alpha_GB=0.5d0
+c alpha_GB=0.01d0
+c alpha_GB1=1.0d0+1.0d0/alpha_GB
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
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
call dyn_ssbond_ene(i,j,evdwij)
c write(iout,*) "PO ZWYKLE", evdwij
+c call flush(iout)
evdw=evdw+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)
+c do k=j+1,iend(i,iint)
+ do k=j+1,nct
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
-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))
+ 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,*) "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
- 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,r_cut_int)
+ 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.eq.0.0d0) cycle
- sssgrad=sscagrad(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+sig0ij
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+c if (energy_dec)
+c & write (iout,*) "rij",1.0d0/rij," rij_shift",rij_shift,
+c & " sig",sig," sig0ij",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
+c if (rij_shift.le.0.0D0) then
+ x0=alpha_GB*(sig-sig0ij)
+ if (energy_dec) write (iout,*) i,j," x0",x0
+ if (rij_shift.le.x0) then
+c sig=2.0d0*sig0ij
+ sigder=-sig*sigsq
+c sigder=0.0d0
+ fac=rij**expon
+ rij12=fac*fac
+c rij12=1.0d0
+ x0=alpha_GB*(sig-sig0ij)
+ facx0=1.0d0/x0**expon
+ facx02=facx0*facx0
+ r012=((1.0d0+alpha_GB)*(sig-sig0ij))**(2*expon)
+ afacx0=aa*facx02
+ bfacx0=bb*facx0
+ abfacx0=afacx0+0.5d0*bfacx0
+ Afac=alpha_GB1*abfacx0
+ Afacsig=0.5d0*alpha_GB1*bfacx0/(sig-sig0ij)
+ BBfac=Afac-(afacx0+bfacx0)
+c BBfac=0.0d0
+ Bfacsig=(-alpha_GB1*(abfacx0+afacx0)+
+ & (afacx0+afacx0+bfacx0))/(sig-sig0ij)
+c Bfacsig=0.0d0
+ Afac=Afac*r012
+ Afacsig=Afacsig*r012
+c Afac=1.0d0
+c Afacsig=0.0d0
+c w(x)=4*eps*((1.0+1.0/alpha_GB)*(y0**12-0.5*y0**6)*(r0/x)**12-(1+1/alpha)*(y0**12-0.5*y0**6)+y0**12-y0**6)
+c eps1=1.0d0
+c eps2rt=1.0d0
+c eps3rt=1.0d0
+ e1 = eps1*eps2rt*eps3rt*Afac*rij12
+ e2 = -eps1*eps2rt*eps3rt*BBfac
+ evdwij = e1+e2
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+c eps2der=0.0d0
+c eps3der=0.0d0
+c eps1_om12=0.0d0
+ evdwij=evdwij*eps2rt*eps3rt
+c Afacsig=0.d0
+c Bfacsig=0.0d0
+ if (lprn) then
+ write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+ write (iout,'(2(a3,i3,2x),18(0pf9.5))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+ & eps1*eps2rt**2*eps3rt**2,om1,om2,om12,
+ & 1.0D0/rij,rij_shift,
+ & evdwij
+ endif
+ if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
+ & 'RE r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
+ evdw=evdw+evdwij*sss
+C Calculate gradient components.
+ e1=e1*eps2rt*eps3rt
+ sigder=-expon*eps1*eps2rt*eps2rt*eps3rt*eps3rt
+ & *(Afacsig*rij12-Bfacsig)*sigder
+ fac=-2.0d0*expon*e1*rij*rij
+c print '(2i4,6f8.4)',i,j,sss,sssgrad*
+c & evdwij,fac,sigma(itypi,itypj),expon
+ fac=fac+evdwij*sssgrad/sss*rij
+c fac=0.0d0
+c write (iout,*) "sigder",sigder," fac",fac," e1",e1,
+c & " e2",e2," sss",sss," sssgrad",sssgrad,"esp123",
+c & eps1*eps2rt**2*eps3rt**2
+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
+C gg_lipi(3)=0.0d0
+C gg_lipj(3)=0.0d0
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+c 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
+c return
+ else
+ rij_shift=1.0D0/rij_shift
+ sigder=-sig*sigsq
c---------------------------------------------------------------
- rij_shift=1.0D0/rij_shift
- fac=rij_shift**expon
+ 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 (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
+ & 'GB r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
+ if (lprn) then
+ write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij
+ 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,'(a,2i5,3f10.5)')
- & 'r sss evdw',i,j,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
-c print '(2i4,6f8.4)',i,j,sss,sssgrad*
-c & evdwij,fac,sigma(itypi,itypj),expon
- fac=fac+evdwij*sssgrad/sss*rij
-c fac=0.0d0
+ 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*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
-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_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
+ endif
C Calculate angular part of the gradient.
c call sc_grad_scale(sss)
- call sc_grad
+ call sc_grad
ENDIF ! dyn_ss
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
C enddo ! zshift
C enddo ! yshift
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.SPLITELE'
- integer xshift,yshift,zshift,subchap
+ double precision boxshift
integer icall
common /srutu/ icall
logical lprn
double precision evdw
- integer itypi,itypj,itypi1,iint,ind
+ 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,dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1
+ & 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
+ gg_lipi=0.0d0
+ gg_lipj=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)
- sss=sscale(1.0d0/rij,r_cut_int)
- if (sss.eq.0.0d0) cycle
- sssgrad=sscagrad(1.0d0/rij,r_cut_int)
+ 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
+ faclip=fac
+ 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+e_augm)*sssgrad/sss*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_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
C Calculate angular part of the gradient.
c call sc_grad_scale(sss)
- call sc_grad
- enddo ! j
- enddo ! iint
+ call sc_grad
+c enddo ! j
+c enddo ! iint
enddo ! i
end
C-----------------------------------------------------------------------------
include 'COMMON.IOUNITS'
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.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),r_cut_int)
sssgrad=sscagrad(sqrt(rij),r_cut_int)
write(iout,*) 'b2=',(b2(k,i-2),k=1,2)
#endif
enddo
- mu=0.0d0
+ mu(:,:nres)=0.0d0
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
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/
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)
+ 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 & (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)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
#ifdef FOURBODY
num_conti=num_cont_hb(i)
#endif
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)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
C xmedi=xmedi+xshift*boxxsize
C ymedi=ymedi+yshift*boxysize
C zmedi=zmedi+zshift*boxzsize
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
& 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 dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi
+ 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/
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
+c write (iout,*) "lipscale",lipscale
iteli=itel(i)
itelj=itel(j)
if (j.eq.i+2 .and. itelj.eq.2) iteli=2
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)
+ 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)
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(dsqrt(rij),r_cut_int)
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-----------------------------------------------------------------------------
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 side-chain vectors.
C
implicit none
+#ifdef MPI
+ include 'mpif.h'
+#endif
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
- integer xshift,yshift,zshift
+ include 'COMMON.TIME1'
double precision ggg(3)
- integer i,iint,j,k,iteli,itypj,subchap
+ 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 xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
- & dist_temp, dist_init
double precision sscale,sscagrad
+ double precision boxshift
+ external boxshift,to_box
+c#ifdef TIMING_ENE
+c double precision time01
+c#endif
evdw2=0.0D0
evdw2_14=0.0d0
c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
C do yshift=-1,1
C do zshift=-1,1
if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
- do i=iatscp_s,iatscp_e
+c do i=iatscp_s,iatscp_e
+ do ikont=g_listscp_start,g_listscp_end
+c#ifdef TIMING_ENE
+c time01=MPI_Wtime()
+c#endif
+ 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"
+!DIR$ INLINE
+ call to_box(xi,yi,zi)
+c do iint=1,nscp_gr(i)
-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)
-
- 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
+!DIR$ INLINE
+ call to_box(xj,yj,zj)
+c#ifdef TIMING_ENE
+c time_escpsetup=time_escpsetup+MPI_Wtime()-time01
+c time01=MPI_Wtime()
+c#endif
+!DIR$ INLINE
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
c print *,xj,yj,zj,'polozenie j'
+c#ifdef TIMING_ENE
+c time_escpsetup=time_escpsetup+MPI_Wtime()-time01
+c time01=MPI_Wtime()
+c#endif
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
c print *,rrij
sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
+c#ifdef TIMING_ENE
+c time_escpcalc=time_escpcalc+MPI_Wtime()-time01
+c#endif
c endif !endif for sscale cutoff
- enddo ! j
+c enddo ! j
- enddo ! iint
+c enddo ! iint
enddo ! i
c enddo !zshift
c enddo !yshift
if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
& iabs(itype(jjj)).eq.1) then
call ssbond_ene(iii,jjj,eij)
- ehpb=ehpb+2*eij
+c ehpb=ehpb+2*eij
+ ehpb=ehpb+eij
endif
cd write (iout,*) "eij",eij
cd & ' waga=',waga,' fac=',fac
enddo
c min_odl=minval(distancek)
- do kk=1,constr_homology
- if(l_homo(kk,ii)) then
- min_odl=distancek(kk)
- exit
- endif
- enddo
- do kk=1,constr_homology
- if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl)
+ if (nexl.gt.0) then
+ min_odl=0.0d0
+ else
+ do kk=1,constr_homology
+ if(l_homo(kk,ii)) then
+ min_odl=distancek(kk)
+ exit
+ endif
+ enddo
+ do kk=1,constr_homology
+ if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl)
& min_odl=distancek(kk)
- enddo
+ enddo
+ endif
c write (iout,* )"min_odl",min_odl
#ifdef DEBUG
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
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.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
+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
+#ifdef DEBUG
+ write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi
+#endif
+ return
+ end
projects / unres.git / blobdiff