include 'COMMON.MD'
include 'COMMON.CONTROL'
include 'COMMON.TIME1'
+ include 'COMMON.SPLITELE'
#ifdef MPI
c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
c & " nfgtasks",nfgtasks
C
C Compute the side-chain and electrostatic interaction energy
C
+C print *,ipot
goto (101,102,103,104,105,106) ipot
C Lennard-Jones potential.
101 call elj(evdw)
goto 107
C Gay-Berne potential (shifted LJ, angular dependence).
104 call egb(evdw)
+C print *,"bylem w egb"
goto 107
C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
105 call egbv(evdw)
eello_turn4=0.0d0
endif
else
-c write (iout,*) "Soft-spheer ELEC potential"
- call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
- & eello_turn4)
+ write (iout,*) "Soft-spheer ELEC potential"
+c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
+c & eello_turn4)
endif
c print *,"Processor",myrank," computed UELEC"
C
C
C Calculate the SC local energy.
C
+C print *,"TU DOCHODZE?"
call esc(escloc)
c print *,"Processor",myrank," computed USC"
C
else
esccor=0.0d0
endif
+C print *,"PRZED MULIt"
c print *,"Processor",myrank," computed Usccorr"
C
C 12/1/95 Multi-body terms
Uconst=0.0d0
Uconst_back=0.0d0
endif
+C 01/27/2015 added by adasko
+C the energy component below is energy transfer into lipid environment
+C based on partition function
+C print *,"przed lipidami"
+ if (wliptran.gt.0) then
+ call Eliptransfer(eliptran)
+ endif
+C print *,"za lipidami"
+ if (AFMlog.gt.0) then
+ call AFMforce(Eafmforce)
+ endif
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
energia(17)=estr
energia(20)=Uconst+Uconst_back
energia(21)=esccor
+ energia(22)=eliptran
+ energia(23)=Eafmforce
c Here are the energies showed per procesor if the are more processors
c per molecule then we sum it up in sum_energy subroutine
c print *," Processor",myrank," calls SUM_ENERGY"
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ eliptran=energia(22)
+ Eafmforce=energia(23)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
+ & +Eafmforce
#endif
energia(0)=etot
c detecting NaNQ
#ifdef MPI
include 'mpif.h'
#endif
- double precision gradbufc(3,maxres),gradbufx(3,maxres),
- & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
+ double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres),
+ & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres)
+ & ,gloc_scbuf(3,-1:maxres)
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
call flush(iout)
#endif
#ifdef SPLITELE
- do i=1,nct
+ do i=0,nct
do j=1,3
gradbufc(j,i)=wsc*gvdwc(j,i)+
& wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
+ & +wliptran*gliptranc(j,i)
+ & +gradafm(j,i)
+
enddo
enddo
#else
- do i=1,nct
+ do i=0,nct
do j=1,3
gradbufc(j,i)=wsc*gvdwc(j,i)+
& wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
+ & +wliptran*gliptranc(j,i)
+ & +gradafm(j,i)
+
enddo
enddo
#endif
enddo
call flush(iout)
#endif
- do i=1,nres
+ do i=0,nres
do j=1,3
gradbufc_sum(j,i)=gradbufc(j,i)
enddo
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,nnt,-1
+ do i=nres-2,-1,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
enddo
call flush(iout)
#endif
- do i=1,nres
+ do i=-1,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,nnt,-1
+ do i=nres-2,-1,-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
+ do i=-1,nct
do j=1,3
#ifdef SPLITELE
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wturn6*gcorr6_turn(j,i)+
& wsccor*gsccorc(j,i)
& +wscloc*gscloc(j,i)
+ & +wliptran*gliptranc(j,i)
+ & +gradafm(j,i)
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& wturn6*gcorr6_turn(j,i)+
& wsccor*gsccorc(j,i)
& +wscloc*gscloc(j,i)
+ & +wliptran*gliptranc(j,i)
+ & +gradafm(j,i)
+
#endif
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wbond*gradbx(j,i)+
& wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
& wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
enddo
#ifdef DEBUG
do i=1,4*nres
glocbuf(i)=gloc(i,icg)
enddo
-#define DEBUG
+c#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc_sc before reduce"
do i=1,nres
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
do i=1,nres
do j=1,3
gloc_scbuf(j,i)=gloc_sc(j,i,icg)
call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
& MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
time_reduce=time_reduce+MPI_Wtime()-time00
-#define DEBUG
+c#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc_sc after reduce"
do i=1,nres
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
#ifdef DEBUG
write (iout,*) "gloc after reduce"
do i=1,4*nres
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ eliptran=energia(22)
+ Eafmforce=energia(23)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
& edihcnstr,ebr*nss,
- & Uconst,etot
+ & Uconst,eliptran,wliptran,Eafmforce,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
+ & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/
& 'ETOT= ',1pE16.6,' (total)')
+
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
& estr,wbond,ebe,wang,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ebr*nss,Uconst,etot
+ & ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
+ & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+C have you changed here?
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
-cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+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
rij=1.0D0/r_inv_ij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+C have you changed here?
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
call sc_angular
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
+C have you changed here?
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
include 'COMMON.SBRIDGE'
logical lprn
+ integer xshift,yshift,zshift
evdw=0.0D0
ccccc energy_dec=.false.
-c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
lprn=.false.
c if (icall.eq.0) lprn=.false.
ind=0
+C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0
+C we have the original box)
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
do i=iatsc_s,iatsc_e
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
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
+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
+
+C xi=xi+xshift*boxxsize
+C yi=yi+yshift*boxysize
+C zi=zi+zshift*boxzsize
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
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)
+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 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)
+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))
+
+c write (iout,'(a7,4f8.3)')
+c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
+ if (sss.gt.0.0d0) then
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
call sc_angular
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+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
+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
+ evdw=evdw+evdwij*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ 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,
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
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
C Calculate angular part of the gradient.
call sc_grad
+ endif
ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
+C enddo ! zshift
+C enddo ! yshift
+C enddo ! xshift
c write (iout,*) "Number of loop steps in EGB:",ind
cccc energy_dec=.false.
return
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
+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
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
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
+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
+C if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5')
+C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
+ 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)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij+e_augm
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ 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),
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
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
include 'COMMON.CALC'
include 'COMMON.IOUNITS'
double precision dcosom1(3),dcosom2(3)
+cc print *,'sss=',sss
eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
enddo
do k=1,3
- gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss
enddo
c write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
- & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
- & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss
c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
cgrad enddo
cgrad enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
enddo
return
end
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
dimension ggg(3)
-cd write(iout,*) 'In EELEC_soft_sphere'
+C write(iout,*) 'In EELEC_soft_sphere'
ees=0.0D0
evdw1=0.0D0
eel_loc=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
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
dxj=dc(1,j)
dyj=dc(2,j)
dzj=dc(3,j)
- xj=c(1,j)+0.5D0*dxj-xmedi
- yj=c(2,j)+0.5D0*dyj-ymedi
- zj=c(3,j)+0.5D0*dzj-zmedi
+ 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
rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(sqrt(rij))
+ sssgrad=sscagrad(sqrt(rij))
if (rij.lt.r0ijsq) then
evdw1ij=0.25d0*(rij-r0ijsq)**2
fac=rij-r0ijsq
evdw1ij=0.0d0
fac=0.0d0
endif
- evdw1=evdw1+evdw1ij
+ evdw1=evdw1+evdw1ij*sss
C
C Calculate contributions to the Cartesian gradient.
C
- ggg(1)=fac*xj
- ggg(2)=fac*yj
- ggg(3)=fac*zj
+ ggg(1)=fac*xj*sssgrad
+ ggg(2)=fac*yj*sssgrad
+ ggg(3)=fac*zj*sssgrad
do k=1,3
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
c write(iout,*) 'b1=',b1(1,i-2)
c write (iout,*) 'theta=', theta(i-1)
enddo
+#else
+ b1(1,i-2)=b(3,iti)
+ b1(2,i-2)=b(5,iti)
+ b2(1,i-2)=b(2,iti)
+ b2(2,i-2)=b(4,iti)
+ b1tilde(1,i-2)=b1(1,i-2)
+ b1tilde(2,i-2)=-b1(2,i-2)
+ b2tilde(1,i-2)=b2(1,i-2)
+ b2tilde(2,i-2)=-b2(2,i-2)
+ EE(1,2,i-2)=eeold(1,2,iti)
+ EE(2,1,i-2)=eeold(2,1,iti)
+ EE(2,2,i-2)=eeold(2,2,iti)
+ EE(1,1,i-2)=eeold(1,1,iti)
+ enddo
+#endif
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do i=3,nres+1
#endif
-#endif
if (i .lt. nres+1) then
sin1=dsin(phi(i))
cos1=dcos(phi(i))
if (i.gt. nnt+2 .and. i.lt.nct+2) then
iti = itortyp(itype(i-2))
else
- iti=ntortyp+1
+ iti=ntortyp
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
iti1 = itortyp(itype(i-1))
else
- iti1=ntortyp+1
+ iti1=ntortyp
endif
cd write (iout,*) '*******i',i,' iti1',iti
cd write (iout,*) 'b1',b1(:,iti)
if (itype(i-1).le.ntyp) then
iti1 = itortyp(itype(i-1))
else
- iti1=ntortyp+1
+ iti1=ntortyp
endif
else
- iti1=ntortyp+1
+ iti1=ntortyp
endif
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
+ include 'COMMON.SPLITELE'
dimension 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),
C
C Loop over i,i+2 and i,i+3 pairs of the peptide groups
C
+C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
do i=iturn3_start,iturn3_end
+ if (i.le.1) cycle
+C write(iout,*) "tu jest i",i
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
- & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+ & .or. itype(i+2).eq.ntyp1
+ & .or. itype(i+3).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ & .or. itype(i+4).eq.ntyp1
+ & ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
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
num_conti=0
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
+ if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+3).eq.ntyp1
- & .or. itype(i+4).eq.ntyp1) cycle
+ & .or. itype(i+4).eq.ntyp1
+ & .or. itype(i+5).eq.ntyp1
+ & .or. itype(i).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ & ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+C Return atom into box, boxxsize is size of box in x dimension
+c 194 continue
+c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
+c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
+C Condition for being inside the proper box
+c if ((xmedi.gt.((0.5d0)*boxxsize)).or.
+c & (xmedi.lt.((-0.5d0)*boxxsize))) then
+c go to 194
+c endif
+c 195 continue
+c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
+c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+C Condition for being inside the proper box
+c if ((ymedi.gt.((0.5d0)*boxysize)).or.
+c & (ymedi.lt.((-0.5d0)*boxysize))) then
+c go to 195
+c endif
+c 196 continue
+c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
+C Condition for being inside the proper box
+c if ((zmedi.gt.((0.5d0)*boxzsize)).or.
+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
+
num_conti=num_cont_hb(i)
c write(iout,*) "JESTEM W PETLI"
call eelecij(i,i+3,ees,evdw1,eel_loc)
& call eturn4(i,eello_turn4)
num_cont_hb(i)=num_conti
enddo ! i
+C Loop over all neighbouring boxes
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
-c do i=7,7
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (i.le.1) cycle
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+ & .or. itype(i+2).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ & ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
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
+C xmedi=xmedi+xshift*boxxsize
+C ymedi=ymedi+yshift*boxysize
+C zmedi=zmedi+zshift*boxzsize
+
+C Return tom into box, boxxsize is size of box in x dimension
+c 164 continue
+c if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
+c if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
+C Condition for being inside the proper box
+c if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
+c & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
+c go to 164
+c endif
+c 165 continue
+c if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
+c if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+C Condition for being inside the proper box
+c if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
+c & (ymedi.lt.((yshift-0.5d0)*boxysize))) then
+c go to 165
+c endif
+c 166 continue
+c if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+c if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
+cC Condition for being inside the proper box
+c if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
+c & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
+c go to 166
+c endif
+
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
-c do j=13,13
-c write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
- if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
+C write (iout,*) i,j
+ if (j.le.1) cycle
+ if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+ & .or.itype(j+2).eq.ntyp1
+ & .or.itype(j-1).eq.ntyp1
+ &) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
enddo ! i
+C enddo ! zshift
+C enddo ! yshift
+C enddo ! xshift
+
c write (iout,*) "Number of loop steps in EELEC:",ind
cd do i=1,nres
cd write (iout,'(i3,3f10.5,5x,3f10.5)')
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
+ include 'COMMON.SPLITELE'
dimension 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),
dx_normj=dc_norm(1,j)
dy_normj=dc_norm(2,j)
dz_normj=dc_norm(3,j)
- xj=c(1,j)+0.5D0*dxj-xmedi
- yj=c(2,j)+0.5D0*dyj-ymedi
- zj=c(3,j)+0.5D0*dzj-zmedi
+C xj=c(1,j)+0.5D0*dxj-xmedi
+C yj=c(2,j)+0.5D0*dyj-ymedi
+C zj=c(3,j)+0.5D0*dzj-zmedi
+ 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
+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))
+c if (sss.gt.0.0d0) then
rrmij=1.0D0/rij
rij=dsqrt(rij)
rmij=1.0D0/rij
ev2=bbb*r6ij
fac3=ael6i*r6ij
fac4=ael3i*r3ij
- evdwij=ev1+ev2
+ evdwij=(ev1+ev2)
el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
el2=fac4*fac
- eesij=el1+el2
+C MARYSIA
+ eesij=(el1+el2)
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
ees=ees+eesij
- evdw1=evdw1+evdwij
+ evdw1=evdw1+evdwij*sss
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,
C Calculate contributions to the Cartesian gradient.
C
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)
+ facvdw=-6*rrmij*(ev1+evdwij)*sss
facel=-3*rrmij*(el1+eesij)
fac1=fac
erij(1)=xj*rmij
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ 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
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
cgrad enddo
cgrad enddo
#else
- facvdw=ev1+evdwij
- facel=el1+eesij
+C MARYSIA
+ facvdw=(ev1+evdwij)*sss
+ facel=(el1+eesij)
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
erij(1)=xj*rmij
cgrad enddo
cgrad enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+ ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
do k=1,3
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
cgrad enddo
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)
+ & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
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)
+ & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+C MARYSIA
+c endif !sscale
IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
& .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
& .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
C Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
& +a33*muij(4)
+c write (iout,*) 'i',i,' j',j,itype(i),itype(j),
+c & ' eel_loc_ij',eel_loc_ij
c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
C Calculate patrial derivative for theta angle
#ifdef NEWCORR
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eelloc',i,j,eel_loc_ij
-c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+c if (eel_loc_ij.ne.0)
+c & write (iout,'(a4,2i4,8f9.5)')'chuj',
+c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
eel_loc=eel_loc+eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
cgrad enddo
C Remaining derivatives of eello
do l=1,3
- gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
- & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
- 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)
- 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)
- 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)
+ 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))
+ 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))
+ 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))
+ 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))
enddo
ENDIF
C Change 12/26/95 to calculate four-body contributions to H-bonding energy
gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
eello_turn4=eello_turn4-(s1+s2+s3)
+c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
+ & 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
+cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
+cd & ' eello_turn4_num',8*eello_turn4_num
#ifdef NEWCORR
gloc(nphi+i,icg)=gloc(nphi+i,icg)
& -(gs13+gsE13+gsEE1)*wturn4
r0_scp=4.5d0
cd print '(a)','Enter ESCP'
cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
do i=iatscp_s,iatscp_e
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))
-
+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 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
+cC 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
+C xi=xi+xshift*boxxsize
+C yi=yi+yshift*boxysize
+C zi=zi+zshift*boxzsize
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
c yj=c(2,nres+j)-yi
c zj=c(3,nres+j)-zi
C Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
+ 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
+C xj=xj-xi
+C yj=yj-yi
+C zj=zj-zi
rij=xj*xj+yj*yj+zj*zj
+
r0ij=r0_scp
r0ijsq=r0ij*r0ij
if (rij.lt.r0ijsq) then
enddo ! iint
enddo ! i
+C enddo !zshift
+C enddo !yshift
+C enddo !xshift
return
end
C-----------------------------------------------------------------------------
include 'COMMON.FFIELD'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
dimension ggg(3)
evdw2=0.0D0
evdw2_14=0.0d0
+c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
cd print '(a)','Enter ESCP'
cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
do i=iatscp_s,iatscp_e
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)
do j=iscpstart(i,iint),iscpend(i,iint)
c yj=c(2,nres+j)-yi
c zj=c(3,nres+j)-zi
C Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
+ 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
+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)))
+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)))
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+e1+e2
+ evdw2_14=evdw2_14+(e1+e2)*sss
endif
evdwij=e1+e2
- evdw2=evdw2+evdwij
+ 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),
& bad(itypj,iteli)
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
C
- fac=-(evdwij+e1)*rrij
+ fac=-(evdwij+e1)*rrij*sss
+ fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- enddo
+c endif !endif for sscale cutoff
+ enddo ! j
enddo ! iint
enddo ! i
+c enddo !zshift
+c enddo !yshift
+c enddo !xshift
do i=1,nct
do j=1,3
gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
- if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
- & iabs(itype(jjj)).eq.1) then
+C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C & iabs(itype(jjj)).eq.1) then
cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
if (.not.dyn_ss .and. i.le.nss) then
C 15/02/13 CC dynamic SSbond - additional check
+ if (ii.gt.nres
+ & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
endif
estr=0.0d0
estr1=0.0d0
do i=ibondp_start,ibondp_end
- if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
- estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
- do j=1,3
- gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
- & *dc(j,i-1)/vbld(i)
- enddo
- if (energy_dec) write(iout,*)
- & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
- 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 gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
+c & *dc(j,i-1)/vbld(i)
+c enddo
+c if (energy_dec) write(iout,*)
+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
+C YES vbldpDUM is the equlibrium length of spring for Dummy atom
+ diff = vbld(i)-vbldpDUM
+ else
+C NO vbldp0 is the equlibrium lenght of spring for peptide group
diff = vbld(i)-vbldp0
- if (energy_dec) write (iout,*)
+ 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
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
enddo
c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
- endif
+c endif
enddo
estr=0.5d0*AKP*estr+estr1
c
nbi=nbondterm(iti)
if (nbi.eq.1) then
diff=vbld(i+nres)-vbldsc0(1,iti)
- if (energy_dec) write (iout,*)
+ if (energy_dec) write (iout,*)
& "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
& AKSC(1,iti),AKSC(1,iti)*diff*diff
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
etheta=0.0D0
c write (*,'(a,i2)') 'EBEND ICG=',icg
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
C Zero the energy function and its derivative at 0 or pi.
call splinthet(theta(i),0.5d0*delta,ss,ssd)
it=itype(i-1)
ichir22=isign(1,itype(i))
endif
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+ if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
y(1)=0.0D0
y(2)=0.0D0
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
z(1)=cos(phii1)
#else
phii1=phi(i+1)
- z(1)=dcos(phii1)
#endif
+ z(1)=dcos(phii1)
z(2)=dsin(phii1)
else
z(1)=0.0D0
bthetk=bthet(k,itype2,ichir21,ichir22)
endif
thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
+c write(iout,*) 'chuj tu', y(k),z(k)
enddo
dthett=thet_pred_mean*ssd
thet_pred_mean=thet_pred_mean*ss+a0thet(it)
& E_theta,E_tc)
endif
etheta=etheta+ethetai
- if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
- & 'ebend',i,ethetai
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
+ & 'ebend',i,ethetai,theta(i),itype(i)
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
C Calculate the contributions to both Gaussian lobes.
C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time)
C The "polynomial part" of the "standard deviation" of this part of
-C the distribution.
+C the distributioni.
+ccc write (iout,*) thetai,thet_pred_mean
sig=polthet(3,it)
do j=2,0,-1
sig=sig*thet_pred_mean+polthet(j,it)
delthe0=thetai-theta0i
term1=-0.5D0*sigcsq*delthec*delthec
term2=-0.5D0*sig0inv*delthe0*delthe0
+C write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0
C Following fuzzy logic is to avoid underflows in dexp and subsequent INFs and
C NaNs in taking the logarithm. We extract the largest exponent which is added
C to the energy (this being the log of the distribution) at the end of energy
C the sum of the contributions from the two lobes and the pre-exponential
C factor. Simple enough, isn't it?
ethetai=(-dlog(termexp)-termm+dlog(termpre))
+C write (iout,*) 'termexp',termexp,termm,termpre,i
C NOW the derivatives!!!
C 6/6/97 Take into account the deformation.
E_theta=(delthec*sigcsq*term1
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
+c print *,i,itype(i-1),itype(i),itype(i-2)
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
+C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
+
if (iabs(itype(i+1)).eq.20) iblock=2
if (iabs(itype(i+1)).ne.20) iblock=1
dethetai=0.0d0
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+ if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
do i=iphi_start,iphi_end
etors_ii=0.0D0
if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
- & .or. itype(i).eq.ntyp1) cycle
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
+ & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+ itori=itortyp(itype(i-2))
+ itori1=itortyp(itype(i-1))
phii=phi(i)
gloci=0.0D0
C Proline-Proline pair is a special case...
c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
- & .or. itype(i).eq.ntyp1) cycle
+C ANY TWO ARE DUMMY ATOMS in row CYCLE
+c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or.
+c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+ if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
+C For introducing the NH3+ and COO- group please check the etor_d for reference
+C and guidance
etors_ii=0.0D0
if (iabs(itype(i)).eq.20) then
iblock=2
etors_d=0.0D0
c write(iout,*) "a tu??"
do i=iphid_start,iphid_end
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
- & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+C ANY TWO ARE DUMMY ATOMS in row CYCLE
+C if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+C & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or.
+C & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1)) .or.
+C & ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle
+ if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
+ & (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
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
itori2=itortyp(itype(i))
gloci2=0.0D0
iblock=1
if (iabs(itype(i+1)).eq.20) iblock=2
+C Iblock=2 Proline type
+C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT
+C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO-
+C if (itype(i+1).eq.ntyp1) iblock=3
+C The problem of NH3+ group can be resolved by adding new parameters please note if there
+C IS or IS NOT need for this
+C IF Yes uncomment below and add to parmread.F appropriate changes and to v1cij and so on
+C is (itype(i-3).eq.ntyp1) ntblock=2
+C ntblock is N-terminal blocking group
C Regular cosine and sine terms
do j=1,ntermd_1(itori,itori1,itori2,iblock)
+C Example of changes for NH3+ blocking group
+C do j=1,ntermd_1(itori,itori1,itori2,iblock,ntblock)
+C v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock)
v1cij=v1c(1,j,itori,itori1,itori2,iblock)
v1sij=v1s(1,j,itori,itori1,itori2,iblock)
v2cij=v1c(2,j,itori,itori1,itori2,iblock)
if (j.lt.nres-1) then
itj1 = itortyp(itype(j+1))
else
- itj1=ntortyp+1
+ itj1=ntortyp
endif
do iii=1,2
dipi(iii,1)=Ub2(iii,i)
if (i.gt.1) then
iti=itortyp(itype(i))
else
- iti=ntortyp+1
+ iti=ntortyp
endif
itk1=itortyp(itype(k+1))
itj=itortyp(itype(j))
if (l.lt.nres-1) then
itl1=itortyp(itype(l+1))
else
- itl1=ntortyp+1
+ itl1=ntortyp
endif
C A1 kernel(j+1) A2T
cd do iii=1,2
if (i.gt.1) then
iti=itortyp(itype(i))
else
- iti=ntortyp+1
+ iti=ntortyp
endif
itk1=itortyp(itype(k+1))
itl=itortyp(itype(l))
if (j.lt.nres-1) then
itj1=itortyp(itype(j+1))
else
- itj1=ntortyp+1
+ itj1=ntortyp
endif
C A2 kernel(j-1)T A1T
call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
if (j.lt.nres-1) then
itj1=itortyp(itype(j+1))
else
- itj1=ntortyp+1
+ itj1=ntortyp
endif
itk=itortyp(itype(k))
itk1=itortyp(itype(k+1))
if (l.lt.nres-1) then
itl1=itortyp(itype(l+1))
else
- itl1=ntortyp+1
+ itl1=ntortyp
endif
#ifdef MOMENT
s1=dip(4,jj,i)*dip(4,kk,k)
if (j.lt.nres-1) then
itj1=itortyp(itype(j+1))
else
- itj1=ntortyp+1
+ itj1=ntortyp
endif
itk=itortyp(itype(k))
if (k.lt.nres-1) then
itk1=itortyp(itype(k+1))
else
- itk1=ntortyp+1
+ itk1=ntortyp
endif
itl=itortyp(itype(l))
if (l.lt.nres-1) then
itl1=itortyp(itype(l+1))
else
- itl1=ntortyp+1
+ itl1=ntortyp
endif
cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l
cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk,
return
end
+CCC----------------------------------------------
+ subroutine Eliptransfer(eliptran)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+C this is done by Adasko
+C print *,"wchodze"
+C structure of box:
+C water
+C--bordliptop-- buffore starts
+C--bufliptop--- here true lipid starts
+C lipid
+C--buflipbot--- lipid ends buffore starts
+C--bordlipbot--buffore ends
+ 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)+c(3,i+1))/2.0d0),boxzsize))
+ 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 ((positi.gt.bordlipbot)
+ &.and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+C print *,"doing sccale for lower part"
+C print *,i,sslip,fracinbuf,ssgradlip
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+C print *, "doing sscalefor top part"
+C print *,i,sslip,fracinbuf,ssgradlip
+ else
+ eliptran=eliptran+pepliptran
+C print *,"I am in true lipid"
+ endif
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ endif
+ enddo
+C print *, "nic nie bylo w lipidzie?"
+C now multiply all by the peptide group transfer factor
+C eliptran=eliptran*pepliptran
+C now the same for side chains
+CV do i=1,1
+ do i=ilip_start,ilip_end
+ if (itype(i).eq.ntyp1) cycle
+ positi=(mod(c(3,i+nres),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+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)
+C print *,positi,bordlipbot,buflipbot
+ if ((positi.gt.bordlipbot)
+ & .and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-
+ &((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i))
+C print *,"I am in true lipid"
+ endif
+ endif ! if in lipid or buffor
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ enddo
+ return
+ end
+C---------------------------------------------------------
+C AFM soubroutine for constant force
+ subroutine AFMforce(Eafmforce)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+ real*8 diffafm(3)
+ dist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ dist=dist+diffafm(i)**2
+ enddo
+ dist=dsqrt(dist)
+ Eafmforce=-forceAFMconst*(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist
+ gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist
+ enddo
+C print *,'AFM',Eafmforce
+ return
+ end
+C AFM soubroutine for constant velocity
+ subroutine AFMvel(Eafmforce)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.NAMES'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.MD'
+ real*8 diffafm(3)
+ dist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ dist=dist+diffafm(i)**2
+ enddo
+ dist=dsqrt(dist)
+ Eafmforce=-(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-(velconst*diffafm(i)/dist-d_t(i,afmend-1))
+ & /d_time
+ gradafm(i,afmbeg-1)=(velconst*diffafm(i)/dist-d_t(i,afmbeg-1))
+ & /d_time
+ enddo
+C print *,'AFM',Eafmforce
+ return
+ end
+
projects / unres.git / blobdiff