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)
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+cmc
+cmc Sep-06: egb takes care of dynamic ss bonds too
+cmc
+c if (dyn_ss) call dyn_set_nss
+
c print *,"Processor",myrank," computed USCSC"
#ifdef TIMING
time01=MPI_Wtime()
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)
+ else if (selfguide.gt.0) then
+ call AFMvel(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"
call sum_energy(energia,.true.)
+ if (dyn_ss) call dyn_set_nss
c print *," Processor",myrank," left SUM_ENERGY"
#ifdef TIMING
time_sumene=time_sumene+MPI_Wtime()-time00
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
#endif
#ifdef MPI
include 'mpif.h'
- double precision gradbufc(3,maxres),gradbufx(3,maxres),
- & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
#endif
+ 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.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)
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+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
yi=c(2,nres+i)
zi=c(3,nres+i)
C Return atom into box, boxxsize is size of box in x dimension
- 134 continue
- if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
- if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+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
- if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 134
- endif
- 135 continue
- if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
- if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+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
- if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
- & (yi.lt.((yshift-0.5d0)*boxysize))) then
- go to 135
- endif
- 136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
+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
- if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 136
+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)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
yj=c(2,nres+j)
zj=c(3,nres+j)
C Return atom J into box the original box
- 137 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+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
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 137
- endif
- 138 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+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
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 138
- endif
- 139 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+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
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 139
+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)
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+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)
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*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=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
- enddo ! zshift
- enddo ! yshift
- enddo ! xshift
+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
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*sss
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ 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*sss
c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
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 Compute the virtual-bond-torsional-angle dependent quantities needed
C to calculate the el-loc multibody terms of various order.
C
+c write(iout,*) 'nphi=',nphi,nres
+#ifdef PARMAT
+ do i=ivec_start+2,ivec_end+2
+#else
+ do i=3,nres+1
+#endif
+#ifdef NEWCORR
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ iti = itortyp(itype(i-2))
+ else
+ iti=ntortyp+1
+ 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
+ endif
+c write(iout,*),i
+ b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0)
+ & +bnew1(2,1,iti)*dsin(theta(i-1))
+ & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0)
+ gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
+ & +bnew1(2,1,iti)*dcos(theta(i-1))
+ & -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
+c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
+c &*(cos(theta(i)/2.0)
+ b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0)
+ & +bnew2(2,1,iti)*dsin(theta(i-1))
+ & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0)
+c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
+c &*(cos(theta(i)/2.0)
+ gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
+ & +bnew2(2,1,iti)*dcos(theta(i-1))
+ & -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
+c if (ggb1(1,i).eq.0.0d0) then
+c write(iout,*) 'i=',i,ggb1(1,i),
+c &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
+c &bnew1(2,1,iti)*cos(theta(i)),
+c &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0
+c endif
+ b1(2,i-2)=bnew1(1,2,iti)
+ gtb1(2,i-2)=0.0
+ b2(2,i-2)=bnew2(1,2,iti)
+ gtb2(2,i-2)=0.0
+ EE(1,1,i-2)=eenew(1,iti)*dcos(theta(i-1))
+ 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)
+ gtEE(1,1,i-2)=-eenew(1,iti)*dsin(theta(i-1))
+ gtEE(1,2,i-2)=0.0d0
+ gtEE(2,2,i-2)=0.0d0
+ gtEE(2,1,i-2)=0.0d0
+c EE(2,2,iti)=0.0d0
+c EE(1,2,iti)=0.5d0*eenew(1,iti)
+c EE(2,1,iti)=0.5d0*eenew(1,iti)
+c b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i))
+c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
+ 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)
+c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
+c write(iout,*) 'b1=',b1(1,i-2)
+c write (iout,*) 'theta=', theta(i-1)
+ enddo
+#else
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ iti = itortyp(itype(i-2))
+ else
+ iti=ntortyp+1
+ 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
+ endif
+ 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
cd write (iout,*) 'Ug',Ug(:,:,i-2)
c if (i .gt. iatel_s+2) then
if (i .gt. nnt+2) then
- call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2))
- call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
+ call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2))
+#ifdef NEWCORR
+ call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2))
+c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
+#endif
+c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti),
+c & EE(1,2,iti),EE(2,2,iti)
+ call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2))
+ call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2))
+c write(iout,*) "Macierz EUG",
+c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2),
+c & eug(2,2,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
& then
call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
enddo
enddo
endif
- call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2))
- call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
+ call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2))
+ call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2))
do k=1,2
muder(k,i-2)=Ub2der(k,i-2)
enddo
iti1=ntortyp
endif
do k=1,2
- mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
+ mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
-cd write (iout,*) 'mu ',mu(:,i-2)
+C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
+c write (iout,*) 'mu ',mu(:,i-2),i-2
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
C Vectors and matrices dependent on a single virtual-bond dihedral.
- call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
+ call matvec2(DD(1,1,iti),b1tilde(1,i-1),auxvec(1))
call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2))
call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2))
call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2))
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),
- & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+ & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij(4)
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
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
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+4).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes by Ana
& .or. itype(i+2).eq.ntyp1
- & .or. itype(i+3).eq.ntyp1
-c & .or. itype(i-1).eq.ntyp1
- & .or. itype(i+4).eq.ntyp1
- & ) cycle
+ & .or. itype(i+3).eq.ntyp1) cycle
+ if(i.gt.1)then
+ if(itype(i-1).eq.ntyp1)cycle
+ end if
+ if(i.LT.nres-3)then
+ if (itype(i+4).eq.ntyp1) cycle
+ end if
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
- 184 continue
- if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
- if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
- if ((xmedi.gt.((0.5d0)*boxxsize)).or.
- & (xmedi.lt.((-0.5d0)*boxxsize))) then
- go to 184
- endif
- 185 continue
- if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
- if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-C Condition for being inside the proper box
- if ((ymedi.gt.((0.5d0)*boxysize)).or.
- & (ymedi.lt.((-0.5d0)*boxysize))) then
- go to 185
- endif
- 186 continue
- if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
- if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-C Condition for being inside the proper box
- if ((zmedi.gt.((0.5d0)*boxzsize)).or.
- & (zmedi.lt.((-0.5d0)*boxzsize))) then
- go to 186
- 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=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
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+5).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes suggested by Ana
& .or. itype(i+3).eq.ntyp1
& .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)
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
- 194 continue
- if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
- if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
+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
- if ((xmedi.gt.((0.5d0)*boxxsize)).or.
- & (xmedi.lt.((-0.5d0)*boxxsize))) then
- go to 194
- endif
- 195 continue
- if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
- if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+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
- if ((ymedi.gt.((0.5d0)*boxysize)).or.
- & (ymedi.lt.((-0.5d0)*boxysize))) then
- go to 195
- endif
- 196 continue
- if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
- if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
+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
- if ((zmedi.gt.((0.5d0)*boxzsize)).or.
- & (zmedi.lt.((-0.5d0)*boxzsize))) then
- go to 196
- endif
+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)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
& call eturn4(i,eello_turn4)
num_cont_hb(i)=num_conti
enddo ! i
C Loop over all neighbouring boxes
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+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
+ if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+2).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes by Ana
& .or. itype(i+2).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
& ) cycle
dxi=dc(1,i)
dyi=dc(2,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
- 164 continue
- if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
- if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
- if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 164
- endif
- 165 continue
- if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
- if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+ 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
- if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
- & (ymedi.lt.((yshift-0.5d0)*boxysize))) then
- go to 165
- endif
- 166 continue
- if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
- if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
+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
- if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 166
- endif
+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 write (iout,*) i,j,itype(i),itype(j)
+C write (iout,*) i,j
+ if (j.le.1) cycle
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((j+2).gt.nres)
+ & .or.((j-1).le.0)
+C end of changes by Ana
& .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
- enddo ! zshift
- enddo ! yshift
- enddo ! xshift
+C enddo ! zshift
+C enddo ! yshift
+C enddo ! xshift
c write (iout,*) "Number of loop steps in EELEC:",ind
cd do i=1,nres
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),
- & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+ & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
+ & gmuij2(4),gmuji2(4)
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
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
- 174 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+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
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 174
- endif
- 175 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+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
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 175
- endif
- 176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+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
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 176
- endif
+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
- xj=xj-xmedi
- yj=yj-ymedi
- zj=zj-zmedi
+C xj=xj-xmedi
+C yj=yj-ymedi
+C zj=zj-zmedi
rij=xj*xj+yj*yj+zj*zj
sss=sscale(sqrt(rij))
C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
C are computed for EVERY pair of non-contiguous peptide groups.
C
+
if (j.lt.nres-1) then
j1=j+1
j2=j-1
j2=j-2
endif
kkk=0
+ lll=0
do k=1,2
do l=1,2
kkk=kkk+1
muij(kkk)=mu(k,i)*mu(l,j)
+c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
+#ifdef NEWCORR
+ gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
+c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
+ gmuij2(kkk)=gUb2(k,i)*mu(l,j)
+ gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
+c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
+ gmuji2(kkk)=mu(k,i)*gUb2(l,j)
+#endif
enddo
enddo
cd write (iout,*) 'EELEC: i',i,' j',j
& +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=',i,j,muij(1),muij(2),muij(3),muij(4)
+C Calculate patrial derivative for theta angle
+#ifdef NEWCORR
+ geel_loc_ij=a22*gmuij1(1)
+ & +a23*gmuij1(2)
+ & +a32*gmuij1(3)
+ & +a33*gmuij1(4)
+c write(iout,*) "derivative over thatai"
+c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
+c & a33*gmuij1(4)
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)+
+ & geel_loc_ij*wel_loc
+c write(iout,*) "derivative over thatai-1"
+c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
+c & a33*gmuij2(4)
+ geel_loc_ij=
+ & a22*gmuij2(1)
+ & +a23*gmuij2(2)
+ & +a32*gmuij2(3)
+ & +a33*gmuij2(4)
+ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
+ & geel_loc_ij*wel_loc
+c Derivative over j residue
+ geel_loc_ji=a22*gmuji1(1)
+ & +a23*gmuji1(2)
+ & +a32*gmuji1(3)
+ & +a33*gmuji1(4)
+c write(iout,*) "derivative over thataj"
+c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+c & a33*gmuji1(4)
+
+ gloc(nphi+j,icg)=gloc(nphi+j,icg)+
+ & geel_loc_ji*wel_loc
+ geel_loc_ji=
+ & +a22*gmuji2(1)
+ & +a23*gmuji2(2)
+ & +a32*gmuji2(3)
+ & +a33*gmuji2(4)
+c write(iout,*) "derivative over thataj-1"
+c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+c & a33*gmuji2(4)
+ gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
+ & geel_loc_ji*wel_loc
+#endif
+cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eelloc',i,j,eel_loc_ij
dimension ggg(3)
double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
& e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
- & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
+ & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),gpizda1(2,2),
+ & gpizda2(2,2),auxgmat1(2,2),auxgmatt1(2,2),
+ & auxgmat2(2,2),auxgmatt2(2,2)
double precision agg(3,4),aggi(3,4),aggi1(3,4),
& aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cd call checkint_turn3(i,a_temp,eello_turn3_num)
call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
+c auxalary matices for theta gradient
+c auxalary matrix for i+1 and constant i+2
+ call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1))
+c auxalary matrix for i+2 and constant i+1
+ call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1))
call transpose2(auxmat(1,1),auxmat1(1,1))
+ call transpose2(auxgmat1(1,1),auxgmatt1(1,1))
+ call transpose2(auxgmat2(1,1),auxgmatt2(1,1))
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
+ call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
+ call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+C Derivatives in theta
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)
+ & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
+ gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
+ & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
+
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
cd write (2,*) 'i,',i,' j',j,'eello_turn3',
dimension ggg(3)
double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
& e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
- & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
+ & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2),
+ & auxgEvec1(2),auxgEvec2(2),auxgEvec3(2),
+ & gte1t(2,2),gte2t(2,2),gte3t(2,2),
+ & gte1a(2,2),gtae3(2,2),gtae3e2(2,2), ae3gte2(2,2),
+ & gtEpizda1(2,2),gtEpizda2(2,2),gtEpizda3(2,2)
double precision agg(3,4),aggi(3,4),aggi1(3,4),
& aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cd call checkint_turn4(i,a_temp,eello_turn4_num)
c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
+c write(iout,*)"WCHODZE W PROGRAM"
a_temp(1,1)=a22
a_temp(1,2)=a23
a_temp(2,1)=a32
call transpose2(EUg(1,1,i+1),e1t(1,1))
call transpose2(Eug(1,1,i+2),e2t(1,1))
call transpose2(Eug(1,1,i+3),e3t(1,1))
+C Ematrix derivative in theta
+ call transpose2(gtEUg(1,1,i+1),gte1t(1,1))
+ call transpose2(gtEug(1,1,i+2),gte2t(1,1))
+ call transpose2(gtEug(1,1,i+3),gte3t(1,1))
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
+c eta1 in derivative theta
+ call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+c auxgvec is derivative of Ub2 so i+3 theta
+ call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1))
+c auxalary matrix of E i+1
+ call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1))
+c s1=0.0
+c gs1=0.0
+ s1=scalar2(b1(1,i+2),auxvec(1))
+c derivative of theta i+2 with constant i+3
+ gs23=scalar2(gtb1(1,i+2),auxvec(1))
+c derivative of theta i+2 with constant i+2
+ gs32=scalar2(b1(1,i+2),auxgvec(1))
+c derivative of E matix in theta of i+1
+ gsE13=scalar2(b1(1,i+2),auxgEvec1(1))
+
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
+c ea31 in derivative theta
+ call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+c auxilary matrix auxgvec of Ub2 with constant E matirx
+ call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
+c auxilary matrix auxgEvec1 of E matix with Ub2 constant
+ call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1))
+
+c s2=0.0
+c gs2=0.0
+ s2=scalar2(b1(1,i+1),auxvec(1))
+c derivative of theta i+1 with constant i+3
+ gs13=scalar2(gtb1(1,i+1),auxvec(1))
+c derivative of theta i+2 with constant i+1
+ gs21=scalar2(b1(1,i+1),auxgvec(1))
+c derivative of theta i+3 with constant i+1
+ gsE31=scalar2(b1(1,i+1),auxgEvec3(1))
+c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),gtb1(1,i+2),
+c & gtb1(1,i+1)
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
+c two derivatives over diffetent matrices
+c gtae3e2 is derivative over i+3
+ call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1))
+c ae3gte2 is derivative over i+2
+ call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
+c three possible derivative over theta E matices
+c i+1
+ call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1))
+c i+2
+ call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1))
+c i+3
+ call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
+
+ gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
+ gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
+ 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
+ gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
+ & -(gs23+gs21+gsEE2)*wturn4
+ gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
+ & -(gs32+gsE31+gsEE3)*wturn4
+c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
+c & gs2
+#endif
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+ & 'eturn4',i,j,-(s1+s2+s3)
+c write (iout,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
+c & ' eello_turn4_num',8*eello_turn4_num
C Derivatives in gamma(i)
call transpose2(EUgder(1,1,i+1),e1tder(1,1))
call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
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)
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))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,1))
call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
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))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3tder(1,1),auxmat(1,1))
call matvec2(auxmat(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1))
call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=agg(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggi(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggi1(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggj(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
a_temp(2,2)=aggj1(l,4)
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+2),auxvec(1))
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
r0_scp=4.5d0
cd print '(a)','Enter ESCP'
cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+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)
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
- 134 continue
- if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
- if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+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
- if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 134
- endif
- 135 continue
- if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
- if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
-C Condition for being inside the proper box
- if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
- & (yi.lt.((yshift-0.5d0)*boxysize))) then
- go to 135
- endif
- 136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
+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
- if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 136
- endif
+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)
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- 174 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 174
- endif
- 175 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+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
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 175
- endif
- 176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+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
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 176
- endif
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+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
enddo ! iint
enddo ! i
- enddo !zshift
- enddo !yshift
- enddo !xshift
+C enddo !zshift
+C enddo !yshift
+C enddo !xshift
return
end
C-----------------------------------------------------------------------------
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
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+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
- 134 continue
- if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
- if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+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
- if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
- & (xi.lt.((xshift-0.5d0)*boxxsize))) then
- go to 134
- endif
- 135 continue
- if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
- if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+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
- if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
- & (yi.lt.((yshift-0.5d0)*boxysize))) then
- go to 135
- endif
- 136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
+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
- if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
- & (zi.lt.((zshift-0.5d0)*boxzsize))) then
- go to 136
- endif
+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)
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- 174 continue
- if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
- if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+ 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
- if ((xj.gt.((0.5d0)*boxxsize)).or.
- & (xj.lt.((-0.5d0)*boxxsize))) then
- go to 174
- endif
- 175 continue
- if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
- if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
- if ((yj.gt.((0.5d0)*boxysize)).or.
- & (yj.lt.((-0.5d0)*boxysize))) then
- go to 175
- endif
- 176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
+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
- if ((zj.gt.((0.5d0)*boxzsize)).or.
- & (zj.lt.((-0.5d0)*boxzsize))) then
- go to 176
- endif
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+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)))
- if (sss.gt.0.0d0) then
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- endif !endif for sscale cutoff
+c endif !endif for sscale cutoff
enddo ! j
enddo ! iint
enddo ! i
- enddo !zshift
- enddo !yshift
- enddo !xshift
+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)
iii=ii
jjj=jj
endif
-cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
+c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+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
cd write (iout,*) "eij",eij
else
C Calculate the distance between the two points and its difference from the
C target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ dd=dist(ii,jj)
+ rdis=dd-dhpb(i)
C Get the force constant corresponding to this distance.
- waga=forcon(i)
+ waga=forcon(i)
C Calculate the contribution to energy.
- ehpb=ehpb+waga*rdis*rdis
+ ehpb=ehpb+waga*rdis*rdis
C
C Evaluate gradient.
C
- fac=waga*rdis/dd
+ fac=waga*rdis/dd
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
- do j=1,3
- ggg(j)=fac*(c(j,jj)-c(j,ii))
- enddo
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
C If this is a SC-SC distance, we need to calculate the contributions to the
C Cartesian gradient in the SC vectors (ghpbx).
- if (iii.lt.ii) then
+ if (iii.lt.ii) then
do j=1,3
ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
- endif
+ endif
cgrad do j=iii,jjj-1
cgrad do k=1,3
cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k)
cgrad enddo
cgrad enddo
- do k=1,3
- ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
- ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
- enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
endif
enddo
ehpb=0.5D0*ehpb
C NO vbldp0 is the equlibrium lenght of spring for peptide group
diff = vbld(i)-vbldp0
endif
- if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+ 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
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=etheta+ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
- gloc(nphi+i-2,icg)=wang*dethetai+ gloc(nphi+i-2,icg)
+ gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg)
enddo
return
end
do iii=1,2
dipi(iii,1)=Ub2(iii,i)
dipderi(iii)=Ub2der(iii,i)
- dipi(iii,2)=b1(iii,iti1)
+ dipi(iii,2)=b1(iii,i+1)
dipj(iii,1)=Ub2(iii,j)
dipderj(iii)=Ub2der(iii,j)
- dipj(iii,2)=b1(iii,itj1)
+ dipj(iii,2)=b1(iii,j+1)
enddo
kkk=0
do iii=1,2
C indluded.
IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN
call transpose2(AEA(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
call transpose2(AEAderg(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
- call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
- call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
+ call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
+ call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
call transpose2(AEA(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,j),AEAb1(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2))
call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2))
call transpose2(AEAderg(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,j),AEAb1derg(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2))
- call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2))
- call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2))
+ call matvec2(AEA(1,1,2),b1(1,l+1),AEAb1(1,2,2))
+ call matvec2(AEAderg(1,1,2),b1(1,l+1),AEAb1derg(1,2,2))
call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2))
call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2))
call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2))
do kkk=1,5
do lll=1,3
call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),
+ call matvec2(auxmat(1,1),b1(1,i),
& AEAb1derx(1,lll,kkk,iii,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),
& AEAb2derx(1,lll,kkk,iii,1,1))
- call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
& AEAb1derx(1,lll,kkk,iii,2,1))
call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
& AEAb2derx(1,lll,kkk,iii,2,1))
call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj),
+ call matvec2(auxmat(1,1),b1(1,j),
& AEAb1derx(1,lll,kkk,iii,1,2))
call matvec2(auxmat(1,1),Ub2(1,j),
& AEAb2derx(1,lll,kkk,iii,1,2))
- call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
& AEAb1derx(1,lll,kkk,iii,2,2))
call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1),
& AEAb2derx(1,lll,kkk,iii,2,2))
IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or.
& (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN
call transpose2(AEA(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
call transpose2(AEAderg(1,1,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
+ call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
- call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
- call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
+ call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
+ call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
call transpose2(AEA(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,j+1),AEAb1(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2))
call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2))
call transpose2(AEAderg(1,1,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2))
+ call matvec2(auxmat(1,1),b1(1,l),AEAb1(1,1,2))
call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2))
- call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2))
- call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2))
+ call matvec2(AEA(1,1,2),b1(1,j+1),AEAb1(1,2,2))
+ call matvec2(AEAderg(1,1,2),b1(1,j+1),AEAb1derg(1,2,2))
call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2))
call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2))
call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2))
do kkk=1,5
do lll=1,3
call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,iti),
+ call matvec2(auxmat(1,1),b1(1,i),
& AEAb1derx(1,lll,kkk,iii,1,1))
call matvec2(auxmat(1,1),Ub2(1,i),
& AEAb2derx(1,lll,kkk,iii,1,1))
- call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
& AEAb1derx(1,lll,kkk,iii,2,1))
call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
& AEAb2derx(1,lll,kkk,iii,2,1))
call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
- call matvec2(auxmat(1,1),b1(1,itl),
+ call matvec2(auxmat(1,1),b1(1,l),
& AEAb1derx(1,lll,kkk,iii,1,2))
call matvec2(auxmat(1,1),Ub2(1,l),
& AEAb2derx(1,lll,kkk,iii,1,2))
- call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1),
+ call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,j+1),
& AEAb1derx(1,lll,kkk,iii,2,2))
call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j),
& AEAb2derx(1,lll,kkk,iii,2,2))
call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
- eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk))
+ eello5_2=scalar2(AEAb1(1,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k))
C Explicit gradient in virtual-dihedral angles.
g_corr5_loc(k-1)=g_corr5_loc(k-1)
vv(2)=pizda(2,1)-pizda(1,2)
if (l.eq.j+1) then
g_corr5_loc(l-1)=g_corr5_loc(l-1)
- & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
+ & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k)))
else
g_corr5_loc(j-1)=g_corr5_loc(j-1)
- & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
+ & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k)))
endif
C Cartesian gradient
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
derx(lll,kkk,iii)=derx(lll,kkk,iii)
- & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk))
+ & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,k))
& -0.5d0*scalar2(vv(1),Ctobr(1,k))
enddo
enddo
call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
- eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl))
+ eello5_4=scalar2(AEAb1(1,2,2),b1(1,l))
& -0.5d0*scalar2(vv(1),Ctobr(1,l))
C Explicit gradient in virtual-dihedral angles.
g_corr5_loc(l-1)=g_corr5_loc(l-1)
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
g_corr5_loc(k-1)=g_corr5_loc(k-1)
- & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl))
+ & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,l))
& -0.5d0*scalar2(vv(1),Ctobr(1,l)))
C Cartesian gradient
do iii=1,2
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
derx(lll,kkk,iii)=derx(lll,kkk,iii)
- & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl))
+ & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,l))
& -0.5d0*scalar2(vv(1),Ctobr(1,l))
enddo
enddo
call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
- eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj))
+ eello5_4=scalar2(AEAb1(1,2,2),b1(1,j))
& -0.5d0*scalar2(vv(1),Ctobr(1,j))
C Explicit gradient in virtual-dihedral angles.
g_corr5_loc(j-1)=g_corr5_loc(j-1)
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
g_corr5_loc(k-1)=g_corr5_loc(k-1)
- & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj))
+ & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,j))
& -0.5d0*scalar2(vv(1),Ctobr(1,j)))
C Cartesian gradient
do iii=1,2
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)
- & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj))
+ & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,j))
& -0.5d0*scalar2(vv(1),Ctobr(1,j))
enddo
enddo
vv1(1)=pizda1(1,1)-pizda1(2,2)
vv1(2)=pizda1(1,2)+pizda1(2,1)
s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
- vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk)
- vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk)
+ vv(1)=AEAb1(1,2,imat)*b1(1,k)-AEAb1(2,2,imat)*b1(2,k)
+ vv(2)=AEAb1(1,2,imat)*b1(2,k)+AEAb1(2,2,imat)*b1(1,k)
s5=scalar2(vv(1),Dtobr2(1,i))
cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5
eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5)
call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1))
vv1(1)=pizda1(1,1)-pizda1(2,2)
vv1(2)=pizda1(1,2)+pizda1(2,1)
- vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk)
- vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk)
+ vv(1)=AEAb1derg(1,2,imat)*b1(1,k)-AEAb1derg(2,2,imat)*b1(2,k)
+ vv(2)=AEAb1derg(1,2,imat)*b1(2,k)+AEAb1derg(2,2,imat)*b1(1,k)
if (l.eq.j+1) then
g_corr6_loc(l-1)=g_corr6_loc(l-1)
& +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i))
vv1(1)=pizda1(1,1)-pizda1(2,2)
vv1(2)=pizda1(1,2)+pizda1(2,1)
s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
- vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk)
- & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk)
- vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk)
- & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk)
+ vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,k)
+ & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,k)
+ vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,k)
+ & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,k)
s5=scalar2(vv(1),Dtobr2(1,i))
derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5)
enddo
C o o C
C \ /l\ /j\ / C
C \ / \ / \ / C
-C o| o | | o |o C
+C o| o | | o |o C
C \ j|/k\| \ |/k\|l C
C \ / \ \ / \ C
C o o C
-C i i C
-C C
+C i i C
+C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
C AL 7/4/01 s1 would occur in the sixth-order moment,
double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C C
+C C
C Parallel Antiparallel C
C C
-C o o C
+C o o C
C /l\ / \ /j\ C
C / \ / \ / \ C
C /| o |o o| o |\ C
#ifdef MOMENT
s1=dip(4,jj,i)*dip(4,kk,k)
#endif
- call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1))
- s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
- call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1))
- s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+ call matvec2(AECA(1,1,1),b1(1,k+1),auxvec(1))
+ s2=0.5d0*scalar2(b1(1,k),auxvec(1))
+ call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1))
+ s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
call transpose2(EE(1,1,itk),auxmat(1,1))
call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
#endif
c eello6_graph3=-s4
C Derivatives in gamma(k-1)
- call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1))
- s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+ call matvec2(AECAderg(1,1,2),b1(1,l+1),auxvec(1))
+ s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k))
g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4)
C Derivatives in gamma(l-1)
- call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1))
- s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
+ call matvec2(AECAderg(1,1,1),b1(1,k+1),auxvec(1))
+ s2=0.5d0*scalar2(b1(1,k),auxvec(1))
call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k)
endif
#endif
- call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+ call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
& auxvec(1))
- s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
- call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
+ s2=0.5d0*scalar2(b1(1,k),auxvec(1))
+ call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
& auxvec(1))
- s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+ s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1),
& pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
& auxvec1(2),auxmat1(2,2)
logical swap
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C C
+C C
C Parallel Antiparallel C
C C
C o o C
C / \ / \ / \ C
C /| o |o o| o |\ C
C \ j|/k\| \ |/k\|l C
-C \ / \ \ / \ C
+C \ / \ \ / \ C
C o \ o \ C
C i i C
-C C
+C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C 4/7/01 AL Component s1 was removed, because it pertains to the respective
call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1))
s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
if (j.eq.l+1) then
- call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+ call matvec2(ADtEA1(1,1,3-imat),b1(1,j+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
else
- call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+ call matvec2(ADtEA1(1,1,3-imat),b1(1,l+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
endif
call transpose2(EUg(1,1,k),auxmat(1,1))
call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1))
#endif
s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1))
if (j.eq.l+1) then
- call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,j+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
else
- call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,l+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
endif
s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i))
if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1))
s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1))
if (j.eq.l+1) then
- call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,j+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
else
- call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+ call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,l+1),auxvec1(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
endif
call transpose2(EUgder(1,1,k),auxmat1(1,1))
call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1))
s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
if (j.eq.l+1) then
call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
- & b1(1,itj1),auxvec(1))
- s3=-0.5d0*scalar2(b1(1,itj),auxvec(1))
+ & b1(1,j+1),auxvec(1))
+ s3=-0.5d0*scalar2(b1(1,j),auxvec(1))
else
call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
- & b1(1,itl1),auxvec(1))
- s3=-0.5d0*scalar2(b1(1,itl),auxvec(1))
+ & b1(1,l+1),auxvec(1))
+ s3=-0.5d0*scalar2(b1(1,l),auxvec(1))
endif
call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1),
& pizda(1,1))
#ifdef MOMENT
call transpose2(AEA(1,1,1),auxmat(1,1))
call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1))
- ss1=scalar2(Ub2(1,i+2),b1(1,itl))
+ ss1=scalar2(Ub2(1,i+2),b1(1,l))
s1 = (auxmat(1,1)+auxmat(2,2))*ss1
#endif
- call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
+ call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1))
- s2 = scalar2(b1(1,itk),vtemp1(1))
+ s2 = scalar2(b1(1,k),vtemp1(1))
#ifdef MOMENT
call transpose2(AEA(1,1,2),atemp(1,1))
call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1))
call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1))
call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1))
call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1))
- ss13 = scalar2(b1(1,itk),vtemp4(1))
+ ss13 = scalar2(b1(1,k),vtemp4(1))
s13 = (gtemp(1,1)+gtemp(2,2))*ss13
#endif
c write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13
#ifdef MOMENT
call transpose2(AEA(1,1,1),auxmatd(1,1))
call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
- ss1d=scalar2(Ub2der(1,i+2),b1(1,itl))
+ ss1d=scalar2(Ub2der(1,i+2),b1(1,l))
s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d
#endif
- call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1))
+ call matvec2(EUgder(1,1,i+2),b1(1,l),vtemp1d(1))
call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1))
- s2d = scalar2(b1(1,itk),vtemp1d(1))
+ s2d = scalar2(b1(1,k),vtemp1d(1))
#ifdef MOMENT
call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1))
s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1))
call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
#endif
- call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1))
+ call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1d(1))
call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1))
- s2d = scalar2(b1(1,itk),vtemp1d(1))
+ s2d = scalar2(b1(1,k),vtemp1d(1))
#ifdef MOMENT
call transpose2(AEA(1,1,2),atempd(1,1))
call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1))
s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
#ifdef MOMENT
call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1))
- ss13d = scalar2(b1(1,itk),vtemp4d(1))
+ ss13d = scalar2(b1(1,k),vtemp4d(1))
s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
#endif
c s1d=0.0d0
call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
#endif
- call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
+ call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1),
& vtemp1d(1))
- s2d = scalar2(b1(1,itk),vtemp1d(1))
+ s2d = scalar2(b1(1,k),vtemp1d(1))
#ifdef MOMENT
call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1))
call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1))
derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d
call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4),
& vtemp4d(1))
- ss13d = scalar2(b1(1,itk),vtemp4d(1))
+ ss13d = scalar2(b1(1,k),vtemp4d(1))
s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d
enddo
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---------------------------------------------------------
+C AFM subroutine with pseudoconstant 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'
+ real*8 diffafm(3)
+C Only for check grad COMMENT if not used for checkgrad
+C totT=3.0d0
+C--------------------------------------------------------
+C print *,"wchodze"
+ 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=0.5d0*forceAFMconst
+ & *(distafminit+totTafm*velAFMconst-dist)**2
+C Eafmforce=-forceAFMconst*(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*
+ &(distafminit+totTafm*velAFMconst-dist)
+ &*diffafm(i)/dist
+ gradafm(i,afmbeg-1)=forceAFMconst*
+ &(distafminit+totTafm*velAFMconst-dist)
+ &*diffafm(i)/dist
+ enddo
+C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
+ return
+ end
projects / unres.git / blobdiff