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
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
logical lprn
integer xshift,yshift,zshift
evdw=0.0D0
go to 135
endif
136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+ if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+ 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 138
endif
139 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+ if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+ 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
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 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)
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(1)=xj*fac
gg(3)=zj*fac
C Calculate angular part of the gradient.
call sc_grad
+ endif
enddo ! j
enddo ! iint
enddo ! i
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)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
- & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss
gvdwx(k,j)=gvdwx(k,j)+gg(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))
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),
go to 185
endif
186 continue
- if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
- if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+ 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 195
endif
196 continue
- if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
- if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+ 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 165
endif
166 continue
- if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
- if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+ if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+ if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
C Condition for being inside the proper box
if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
& (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
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),
go to 175
endif
176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+ if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+ 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
yj=yj-ymedi
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 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)
+ fac=-3*rrmij*(facvdw+facvdw+facel)+sssgrad*rmij*evdwij
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*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*sss
+ ggg(2)=facvdw*yj*sss
+ ggg(3)=facvdw*zj*sss
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
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
& 'eelloc',i,j,eel_loc_ij
- if (eel_loc_ij.ne.0)
- & write (iout,'(a4,2i4,8f9.5)')'chuj',
- & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
+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
go to 135
endif
136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+ if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+ 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 175
endif
176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+ if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+ 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
yj=yj-yi
zj=zj-zi
rij=xj*xj+yj*yj+zj*zj
+
r0ij=r0_scp
r0ijsq=r0ij*r0ij
if (rij.lt.r0ijsq) then
include 'COMMON.FFIELD'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
dimension ggg(3)
evdw2=0.0D0
evdw2_14=0.0d0
go to 135
endif
136 continue
- if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
- if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+ if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+ 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 175
endif
176 continue
- if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
- if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+ if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+ 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
yj=yj-yi
zj=zj-zi
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ sss=sscale(1.0d0/(dsqrt(rrij)))
+ 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)
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
+ endif !endif for sscale cutoff
+ enddo ! j
enddo ! iint
enddo ! i
etheta=0.0D0
c write (*,'(a,i2)') 'EBEND ICG=',icg
do i=ithet_start,ithet_end
+ print *,i,itype(i-1),itype(i),itype(i-2)
if (itype(i-1).eq.ntyp1) cycle
+ print *,'wchodze',itype(i-1)
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)
& '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)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
enddo
C Ufff.... We've done all this!!!
return
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 distributioni.
- write (iout,*) thetai,thet_pred_mean
+ccc write (iout,*) thetai,thet_pred_mean
sig=polthet(3,it)
do j=2,0,-1
sig=sig*thet_pred_mean+polthet(j,it)
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
+c print *,i,itype(i-1),itype(i),itype(i-2)
if ((itype(i-1).eq.ntyp1)) cycle
if (iabs(itype(i+1)).eq.20) iblock=2
if (iabs(itype(i+1)).ne.20) iblock=1
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)=wang*dethetai+ gloc(nphi+i-2,icg)
enddo
return
end
etors=0.0D0
do i=iphi_start,iphi_end
C ANY TWO ARE DUMMY ATOMS in row CYCLE
- if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
- & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or.
- & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) 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-3).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+ & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
C For introducing the NH3+ and COO- group please check the etor_d for reference
C and guidance
etors_ii=0.0D0
projects / unres.git / blobdiff