+! Check the gradient of the virtual-bond and SC vectors in the internal
+! coordinates.
+!
+ aincr=1.0d-6
+ aincr2=5.0d-7
+ call cartder
+ write (iout,'(a)') '**************** dx/dalpha'
+ write (iout,'(a)')
+ do i=2,nres-1
+ alphi=alph(i)
+ alph(i)=alph(i)+aincr
+ do k=1,3
+ temp(k,i)=dc(k,nres+i)
+ enddo
+ call chainbuild
+ do k=1,3
+ gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
+ xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr))
+ enddo
+ write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') &
+ i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3)
+ write (iout,'(a)')
+ alph(i)=alphi
+ call chainbuild
+ enddo
+ write (iout,'(a)')
+ write (iout,'(a)') '**************** dx/domega'
+ write (iout,'(a)')
+ do i=2,nres-1
+ omegi=omeg(i)
+ omeg(i)=omeg(i)+aincr
+ do k=1,3
+ temp(k,i)=dc(k,nres+i)
+ enddo
+ call chainbuild
+ do k=1,3
+ gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
+ xx(k)=dabs((gg(k)-dxds(k+3,i))/ &
+ (aincr*dabs(dxds(k+3,i))+aincr))
+ enddo
+ write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') &
+ i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3)
+ write (iout,'(a)')
+ omeg(i)=omegi
+ call chainbuild
+ enddo
+ write (iout,'(a)')
+ write (iout,'(a)') '**************** dx/dtheta'
+ write (iout,'(a)')
+ do i=3,nres
+ theti=theta(i)
+ theta(i)=theta(i)+aincr
+ do j=i-1,nres-1
+ do k=1,3
+ temp(k,j)=dc(k,nres+j)
+ enddo
+ enddo
+ call chainbuild
+ do j=i-1,nres-1
+ ii = indmat(i-2,j)
+! print *,'i=',i-2,' j=',j-1,' ii=',ii
+ do k=1,3
+ gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
+ xx(k)=dabs((gg(k)-dxdv(k,ii))/ &
+ (aincr*dabs(dxdv(k,ii))+aincr))
+ enddo
+ write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') &
+ i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3)
+ write(iout,'(a)')
+ enddo
+ write (iout,'(a)')
+ theta(i)=theti
+ call chainbuild
+ enddo
+ write (iout,'(a)') '***************** dx/dphi'
+ write (iout,'(a)')
+ do i=4,nres
+ phi(i)=phi(i)+aincr
+ do j=i-1,nres-1
+ do k=1,3
+ temp(k,j)=dc(k,nres+j)
+ enddo
+ enddo
+ call chainbuild
+ do j=i-1,nres-1
+ ii = indmat(i-2,j)
+! print *,'ii=',ii
+ do k=1,3
+ gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
+ xx(k)=dabs((gg(k)-dxdv(k+3,ii))/ &
+ (aincr*dabs(dxdv(k+3,ii))+aincr))
+ enddo
+ write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') &
+ i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3)
+ write(iout,'(a)')
+ enddo
+ phi(i)=phi(i)-aincr
+ call chainbuild
+ enddo
+ write (iout,'(a)') '****************** ddc/dtheta'
+ do i=1,nres-2
+ thet=theta(i+2)
+ theta(i+2)=thet+aincr
+ do j=i,nres
+ do k=1,3
+ temp(k,j)=dc(k,j)
+ enddo
+ enddo
+ call chainbuild
+ do j=i+1,nres-1
+ ii = indmat(i,j)
+! print *,'ii=',ii
+ do k=1,3
+ gg(k)=(dc(k,j)-temp(k,j))/aincr
+ xx(k)=dabs((gg(k)-dcdv(k,ii))/ &
+ (aincr*dabs(dcdv(k,ii))+aincr))
+ enddo
+ write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') &
+ i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3)
+ write (iout,'(a)')
+ enddo
+ do j=1,nres
+ do k=1,3
+ dc(k,j)=temp(k,j)
+ enddo
+ enddo
+ theta(i+2)=thet
+ enddo
+ write (iout,'(a)') '******************* ddc/dphi'
+ do i=1,nres-3
+ phii=phi(i+3)
+ phi(i+3)=phii+aincr
+ do j=1,nres
+ do k=1,3
+ temp(k,j)=dc(k,j)
+ enddo
+ enddo
+ call chainbuild
+ do j=i+2,nres-1
+ ii = indmat(i+1,j)
+! print *,'ii=',ii
+ do k=1,3
+ gg(k)=(dc(k,j)-temp(k,j))/aincr
+ xx(k)=dabs((gg(k)-dcdv(k+3,ii))/ &
+ (aincr*dabs(dcdv(k+3,ii))+aincr))
+ enddo
+ write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') &
+ i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3)
+ write (iout,'(a)')
+ enddo
+ do j=1,nres
+ do k=1,3
+ dc(k,j)=temp(k,j)
+ enddo
+ enddo
+ phi(i+3)=phii
+ enddo
+ return
+ end subroutine check_cartgrad
+!-----------------------------------------------------------------------------
+ subroutine check_ecart
+! Check the gradient of the energy in Cartesian coordinates.
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.VAR'
+! include 'COMMON.CONTACTS'
+ use comm_srutu
+!#ifdef LBFGS
+! use minimm, only: funcgrad
+!#endif
+!el integer :: icall
+!el common /srutu/ icall
+! real(kind=8) :: funcgrad
+ real(kind=8),dimension(6) :: ggg
+ real(kind=8),dimension(3) :: cc,xx,ddc,ddx
+ real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(6,nres) :: grad_s
+ real(kind=8),dimension(0:n_ene) :: energia,energia1
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+!EL external fdum
+ integer :: nf,i,j,k
+ real(kind=8) :: aincr,etot,etot1,ff
+ icg=1
+ nf=0
+ nfl=0
+ call zerograd
+ aincr=1.0D-5
+ print '(a)','CG processor',me,' calling CHECK_CART.',aincr
+ nf=0
+ icall=0
+ call geom_to_var(nvar,x)
+ call etotal(energia)
+ etot=energia(0)
+#ifdef LBFGS
+ ff=funcgrad(x,g)
+#else
+!el call enerprint(energia)
+ call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
+#endif
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gradc(j,i,icg)
+ grad_s(j+3,i)=gradx(j,i,icg)
+ enddo
+ enddo
+ call flush(iout)
+ write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
+ do i=1,nres
+ do j=1,3
+ xx(j)=c(j,i+nres)
+ ddc(j)=dc(j,i)
+ ddx(j)=dc(j,i+nres)
+ enddo
+ do j=1,3
+ dc(j,i)=dc(j,i)+aincr
+ do k=i+1,nres
+ c(j,k)=c(j,k)+aincr
+ c(j,k+nres)=c(j,k+nres)+aincr
+ enddo
+ call zerograd
+ call etotal(energia1)
+ etot1=energia1(0)
+ ggg(j)=(etot1-etot)/aincr
+ dc(j,i)=ddc(j)
+ do k=i+1,nres
+ c(j,k)=c(j,k)-aincr
+ c(j,k+nres)=c(j,k+nres)-aincr
+ enddo
+ enddo
+ do j=1,3
+ c(j,i+nres)=c(j,i+nres)+aincr
+ dc(j,i+nres)=dc(j,i+nres)+aincr
+ call zerograd
+ call etotal(energia1)
+ etot1=energia1(0)
+ ggg(j+3)=(etot1-etot)/aincr
+ c(j,i+nres)=xx(j)
+ dc(j,i+nres)=ddx(j)
+ enddo
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
+ enddo
+ return
+ end subroutine check_ecart
+#ifdef CARGRAD
+!-----------------------------------------------------------------------------
+ subroutine check_ecartint
+! Check the gradient of the energy in Cartesian coordinates.
+ use io_base, only: intout
+ use MD_data, only: iset
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.VAR'
+! include 'COMMON.CONTACTS'
+! include 'COMMON.MD'
+! include 'COMMON.LOCAL'
+! include 'COMMON.SPLITELE'
+ use comm_srutu
+!el integer :: icall
+!el common /srutu/ icall
+ real(kind=8),dimension(6) :: ggg,ggg1
+ real(kind=8),dimension(3) :: cc,xx,ddc,ddx,ddc1,ddcn
+ real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(3) :: dcnorm_safe1,dcnorm_safe2,dxnorm_safe
+ real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
+ real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
+ real(kind=8),dimension(0:n_ene) :: energia,energia1
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+!EL external fdum
+ integer :: i,j,k,nf
+ real(kind=8) :: rlambd,aincr,etot,etot1,etot11,etot12,etot2,&
+ etot21,etot22
+ r_cut=2.0d0
+ rlambd=0.3d0
+ icg=1
+ nf=0
+ nfl=0
+ if (iset.eq.0) iset=1
+ call intout
+! call intcartderiv
+! call checkintcartgrad
+ call zerograd
+ aincr=graddelta
+ write(iout,*) 'Calling CHECK_ECARTINT.,kupa'
+ nf=0
+ icall=0
+ call geom_to_var(nvar,x)
+ write (iout,*) "split_ene ",split_ene
+ call flush(iout)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia)
+ etot=energia(0)
+ call cartgrad
+#ifdef FIVEDIAG
+ call grad_transform
+#endif
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ write(iout,*) "before movement analytical gradient"
+
+ enddo
+ enddo
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+
+ else
+!- split gradient check
+ call zerograd
+ call etotal_long(energia)
+!el call enerprint(energia)
+ call cartgrad
+#ifdef FIVEDIAG
+ call grad_transform
+#endif
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ call zerograd
+ call etotal_short(energia)
+ call enerprint(energia)
+ call cartgrad
+#ifdef FIVEDIAG
+ call grad_transform
+#endif
+
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s1(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s1(j,i)=gcart(j,i)
+ grad_s1(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ endif
+ write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
+#ifdef FIVEDIAG
+ do i=1,nres
+#else
+ do i=nnt,nct
+#endif
+ do j=1,3
+ if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
+ if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
+ ddc(j)=c(j,i)
+ ddx(j)=c(j,i+nres)
+ dcnorm_safe1(j)=dc_norm(j,i-1)
+ dcnorm_safe2(j)=dc_norm(j,i)
+ dxnorm_safe(j)=dc_norm(j,i+nres)
+ enddo
+ do j=1,3
+ c(j,i)=ddc(j)+aincr
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+ etot1=energia1(0)
+! write (iout,*) "ij",i,j," etot1",etot1
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+ endif
+!- end split gradient
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
+ c(j,i)=ddc(j)-aincr
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+ etot2=energia1(0)
+! write (iout,*) "ij",i,j," etot2",etot2
+ ggg(j)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+! write (iout,*) "etot21",etot21," etot22",etot22
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ c(j,i)=ddc(j)
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i-1)=dcnorm_safe1(j)
+ dc_norm(j,i)=dcnorm_safe2(j)
+ dc_norm(j,i+nres)=dxnorm_safe(j)
+ enddo
+ do j=1,3
+ c(j,i+nres)=ddx(j)+aincr
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+ etot1=energia1(0)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+ endif
+!- end split gradient
+ c(j,i+nres)=ddx(j)-aincr
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+ etot2=energia1(0)
+ ggg(j+3)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j+3)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j+3)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ c(j,i+nres)=ddx(j)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dxnorm_safe(j)
+ call int_from_cart1(.false.)
+ enddo
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
+ if (split_ene) then
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),&
+ k=1,6)
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),&
+ ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
+ endif
+ enddo
+ return
+ end subroutine check_ecartint
+#else
+!-----------------------------------------------------------------------------
+ subroutine check_ecartint
+! Check the gradient of the energy in Cartesian coordinates.
+ use io_base, only: intout
+ use MD_data, only: iset
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.VAR'
+! include 'COMMON.CONTACTS'
+! include 'COMMON.MD'
+! include 'COMMON.LOCAL'
+! include 'COMMON.SPLITELE'
+ use comm_srutu
+!el integer :: icall
+!el common /srutu/ icall
+ real(kind=8),dimension(6) :: ggg,ggg1
+ real(kind=8),dimension(3) :: cc,xx,ddc,ddx
+ real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(3) :: dcnorm_safe,dxnorm_safe
+ real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
+ real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
+ real(kind=8),dimension(0:n_ene) :: energia,energia1
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+!EL external fdum
+ integer :: i,j,k,nf
+ real(kind=8) :: rlambd,aincr,etot,etot1,etot11,etot12,etot2,&
+ etot21,etot22
+ r_cut=2.0d0
+ rlambd=0.3d0
+ icg=1
+ nf=0
+ nfl=0
+ if (iset.eq.0) iset=1
+ call intout
+! call intcartderiv
+! call checkintcartgrad
+ call zerograd
+ aincr=1.0D-6
+ write(iout,*) 'Calling CHECK_ECARTINT.',aincr
+ nf=0
+ icall=0
+ call geom_to_var(nvar,x)
+ if (.not.split_ene) then
+ call etotal(energia)
+ etot=energia(0)
+! call enerprint(energia)
+ call cartgrad
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ grad_s(j+3,0)=gxcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ write(iout,*) "before movement analytical gradient"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+
+ else
+!- split gradient check
+ call zerograd
+ call etotal_long(energia)
+!el call enerprint(energia)
+ call cartgrad
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ call zerograd
+ call etotal_short(energia)
+!el call enerprint(energia)
+ call cartgrad
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s1(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s1(j,i)=gcart(j,i)
+ grad_s1(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ endif
+ write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
+ do i=0,nres
+ do j=1,3
+ xx(j)=c(j,i+nres)
+ ddc(j)=dc(j,i)
+ ddx(j)=dc(j,i+nres)
+ do k=1,3
+ dcnorm_safe(k)=dc_norm(k,i)
+ dxnorm_safe(k)=dc_norm(k,i+nres)
+ enddo
+ enddo
+ do j=1,3
+ dc(j,i)=ddc(j)+aincr
+ call chainbuild_cart
+#ifdef MPI
+! Broadcast the order to compute internal coordinates to the slaves.
+! if (nfgtasks.gt.1)
+! & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
+#endif
+! call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+ etot1=energia1(0)
+! call enerprint(energia1)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+! write (iout,*) "etot11",etot11," etot12",etot12
+ endif
+!- end split gradient
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
+ dc(j,i)=ddc(j)-aincr
+ call chainbuild_cart
+! call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+! call enerprint(energia1)
+ etot2=energia1(0)
+ ggg(j)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+! write (iout,*) "etot21",etot21," etot22",etot22
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ dc(j,i)=ddc(j)
+ call chainbuild_cart
+ enddo
+ do j=1,3
+ dc(j,i+nres)=ddx(j)+aincr
+ call chainbuild_cart
+! write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
+! write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
+! write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
+! write (iout,*) "dxnormnorm",dsqrt(
+! & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
+! write (iout,*) "dxnormnormsafe",dsqrt(
+! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
+! write (iout,*)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+! call enerprint(energia1)
+ etot1=energia1(0)
+! print *,"ene",energia1(0),energia1(57)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+ endif
+!- end split gradient
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
+ dc(j,i+nres)=ddx(j)-aincr
+ call chainbuild_cart
+! write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
+! write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
+! write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
+! write (iout,*)
+! write (iout,*) "dxnormnorm",dsqrt(
+! & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
+! write (iout,*) "dxnormnormsafe",dsqrt(
+! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
+ if (.not.split_ene) then
+ call zerograd
+ call etotal(energia1)
+ etot2=energia1(0)
+! call enerprint(energia1)
+! print *,"ene",energia1(0),energia1(57)
+ ggg(j+3)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j+3)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j+3)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ dc(j,i+nres)=ddx(j)
+ call chainbuild_cart
+ enddo
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
+ if (split_ene) then
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),&
+ k=1,6)
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),&
+ ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
+ endif
+ enddo
+ return
+ end subroutine check_ecartint
+#endif
+!-----------------------------------------------------------------------------
+ subroutine check_eint
+! Check the gradient of energy in internal coordinates.
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.VAR'
+! include 'COMMON.GEO'
+ use comm_srutu
+!#ifdef LBFGS
+! use minimm, only : funcgrad
+!#endif
+!el integer :: icall
+!el common /srutu/ icall
+! real(kind=8) :: funcgrad
+ real(kind=8),dimension(6*nres) :: x,gana,gg !(maxvar) (maxvar=6*maxres)
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+ real(kind=8),dimension(0:n_ene) :: energia,energia1,energia2
+ character(len=6) :: key
+!EL external fdum
+ integer :: i,ii,nf
+ real(kind=8) :: xi,aincr,etot,etot1,etot2,ff
+ call zerograd
+ aincr=1.0D-7
+ print '(a)','Calling CHECK_INT.'
+ nf=0
+ nfl=0
+ icg=1
+ call geom_to_var(nvar,x)
+ call var_to_geom(nvar,x)
+ call chainbuild
+ icall=1
+! print *,'ICG=',ICG
+ call etotal(energia)
+ etot = energia(0)
+!el call enerprint(energia)
+! print *,'ICG=',ICG
+#ifdef MPL
+ if (MyID.ne.BossID) then
+ call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
+ nf=x(nvar+1)
+ nfl=x(nvar+2)
+ icg=x(nvar+3)
+ endif
+#endif
+ nf=1
+ nfl=3
+#ifdef LBFGS
+ ff=funcgrad(x,gana)
+#else
+
+!d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
+ call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
+!d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp
+#endif
+ icall=1
+ do i=1,nvar
+ xi=x(i)
+ x(i)=xi-0.5D0*aincr
+ call var_to_geom(nvar,x)
+ call chainbuild
+ call etotal(energia1)
+ etot1=energia1(0)
+ x(i)=xi+0.5D0*aincr
+ call var_to_geom(nvar,x)
+ call chainbuild
+ call etotal(energia2)
+ etot2=energia2(0)
+ gg(i)=(etot2-etot1)/aincr
+ write (iout,*) i,etot1,etot2
+ x(i)=xi
+ enddo
+ write (iout,'(/2a)')' Variable Numerical Analytical',&
+ ' RelDiff*100% '
+ do i=1,nvar
+ if (i.le.nphi) then
+ ii=i
+ key = ' phi'
+ else if (i.le.nphi+ntheta) then
+ ii=i-nphi
+ key=' theta'
+ else if (i.le.nphi+ntheta+nside) then
+ ii=i-(nphi+ntheta)
+ key=' alpha'
+ else
+ ii=i-(nphi+ntheta+nside)
+ key=' omega'
+ endif
+ write (iout,'(i3,a,i3,3(1pd16.6))') &
+ i,key,ii,gg(i),gana(i),&
+ 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
+ enddo
+ return
+ end subroutine check_eint
+!-----------------------------------------------------------------------------
+! econstr_local.F
+!-----------------------------------------------------------------------------
+ subroutine Econstr_back
+! MD with umbrella_sampling using Wolyne's distance measure as a constraint
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.VAR'
+! include 'COMMON.MD'
+ use MD_data
+!#ifndef LANG0
+! include 'COMMON.LANGEVIN'
+!#else
+! include 'COMMON.LANGEVIN.lang0'
+!#endif
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+! include 'COMMON.TIME1'
+ integer :: i,j,ii,k
+ real(kind=8) :: utheta_i,dtheta_i,ugamma_i,dgamma_i,dxx,dyy,dzz
+
+ if(.not.allocated(utheta)) allocate(utheta(nfrag_back))
+ if(.not.allocated(ugamma)) allocate(ugamma(nfrag_back))
+ if(.not.allocated(uscdiff)) allocate(uscdiff(nfrag_back))
+
+ Uconst_back=0.0d0
+ do i=1,nres
+ dutheta(i)=0.0d0
+ dugamma(i)=0.0d0
+ do j=1,3
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
+ enddo
+ do i=1,nfrag_back
+ ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset)
+!
+! Deviations from theta angles
+!
+ utheta_i=0.0d0
+ do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset)
+ dtheta_i=theta(j)-thetaref(j)
+ utheta_i=utheta_i+0.5d0*dtheta_i*dtheta_i
+ dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
+ enddo
+ utheta(i)=utheta_i/(ii-1)
+!
+! Deviations from gamma angles
+!
+ ugamma_i=0.0d0
+ do j=ifrag_back(1,i,iset)+3,ifrag_back(2,i,iset)
+ dgamma_i=pinorm(phi(j)-phiref(j))
+! write (iout,*) j,phi(j),phi(j)-phiref(j)
+ ugamma_i=ugamma_i+0.5d0*dgamma_i*dgamma_i
+ dugamma(j-3)=dugamma(j-3)+wfrag_back(2,i,iset)*dgamma_i/(ii-2)
+! write (iout,*) i,j,dgamma_i,wfrag_back(2,i,iset),dugamma(j-3)
+ enddo
+ ugamma(i)=ugamma_i/(ii-2)
+!
+! Deviations from local SC geometry
+!
+ uscdiff(i)=0.0d0
+ do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1
+ dxx=xxtab(j)-xxref(j)
+ dyy=yytab(j)-yyref(j)
+ dzz=zztab(j)-zzref(j)
+ uscdiff(i)=uscdiff(i)+dxx*dxx+dyy*dyy+dzz*dzz
+ do k=1,3
+ duscdiff(k,j-1)=duscdiff(k,j-1)+wfrag_back(3,i,iset)* &
+ (dXX_C1tab(k,j)*dxx+dYY_C1tab(k,j)*dyy+dZZ_C1tab(k,j)*dzz)/ &
+ (ii-1)
+ duscdiff(k,j)=duscdiff(k,j)+wfrag_back(3,i,iset)* &
+ (dXX_Ctab(k,j)*dxx+dYY_Ctab(k,j)*dyy+dZZ_Ctab(k,j)*dzz)/ &
+ (ii-1)
+ duscdiffx(k,j)=duscdiffx(k,j)+wfrag_back(3,i,iset)* &
+ (dXX_XYZtab(k,j)*dxx+dYY_XYZtab(k,j)*dyy+dZZ_XYZtab(k,j)*dzz) &
+ /(ii-1)
+ enddo
+! write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
+! & xxref(j),yyref(j),zzref(j)
+ enddo
+ uscdiff(i)=0.5d0*uscdiff(i)/(ii-1)
+! write (iout,*) i," uscdiff",uscdiff(i)
+!
+! Put together deviations from local geometry
+!
+ Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+ &
+ wfrag_back(2,i,iset)*ugamma(i)+wfrag_back(3,i,iset)*uscdiff(i)
+! write(iout,*) "i",i," utheta",utheta(i)," ugamma",ugamma(i),
+! & " uconst_back",uconst_back
+ utheta(i)=dsqrt(utheta(i))
+ ugamma(i)=dsqrt(ugamma(i))
+ uscdiff(i)=dsqrt(uscdiff(i))
+ enddo
+ return
+ end subroutine Econstr_back
+!-----------------------------------------------------------------------------
+! energy_p_new-sep_barrier.F
+!-----------------------------------------------------------------------------
+ real(kind=8) function sscale(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscale=1.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale=0d0
+ endif
+ return
+ end function sscale
+ real(kind=8) function sscale_grad(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscale_grad=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscale_grad=0d0
+ endif
+ return
+ end function sscale_grad
+!SCALINING MARTINI
+ real(kind=8) function sscale_martini(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+! print *,"here2",r_cut_mart,r
+ if(r.lt.r_cut_mart-rlamb_mart) then
+ sscale_martini=1.0d0
+ else if(r.le.r_cut_mart.and.r.ge.r_cut_mart-rlamb_mart) then
+ gamm=(r-(r_cut_mart-rlamb_mart))/rlamb_mart
+ sscale_martini=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale_martini=0.0d0
+ endif
+ return
+ end function sscale_martini
+ real(kind=8) function sscale_grad_martini(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut_mart-rlamb_mart) then
+ sscale_grad_martini=0.0d0
+ else if(r.le.r_cut_mart.and.r.ge.r_cut_mart-rlamb_mart) then
+ gamm=(r-(r_cut_mart-rlamb_mart))/rlamb_mart
+ sscale_grad_martini=gamm*(6*gamm-6.0d0)/rlamb_mart
+ else
+ sscale_grad_martini=0.0d0
+ endif
+ return
+ end function sscale_grad_martini
+ real(kind=8) function sscale_martini_angle(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm,r_cut_angle,rlamb_angle
+! print *,"here2",r_cut_angle,r
+ r_cut_angle=3.12d0
+ rlamb_angle=0.1d0
+ if(r.lt.r_cut_angle-rlamb_angle) then
+ sscale_martini_angle=1.0d0
+ else if(r.le.r_cut_angle.and.r.ge.r_cut_angle-rlamb_angle) then
+ gamm=(r-(r_cut_angle-rlamb_angle))/rlamb_angle
+ sscale_martini_angle=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale_martini_angle=0.0d0
+ endif
+ return
+ end function sscale_martini_angle
+ real(kind=8) function sscale_grad_martini_angle(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm,r_cut_angle,rlamb_angle
+ r_cut_angle=3.12d0
+ rlamb_angle=0.1d0
+ if(r.lt.r_cut_angle-rlamb_angle) then
+ sscale_grad_martini_angle=0.0d0
+ else if(r.le.r_cut_angle.and.r.ge.r_cut_angle-rlamb_angle) then
+ gamm=(r-(r_cut_angle-rlamb_angle))/rlamb_angle
+ sscale_grad_martini_angle=gamm*(6*gamm-6.0d0)/rlamb_angle
+ else
+ sscale_grad_martini_angle=0.0d0
+ endif
+ return
+ end function sscale_grad_martini_angle
+
+
+!!!!!!!!!! PBCSCALE
+ real(kind=8) function sscale_ele(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscale_ele=1.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscale_ele=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale_ele=0d0
+ endif
+ return
+ end function sscale_ele
+
+ real(kind=8) function sscagrad_ele(r)
+ real(kind=8) :: r,gamm
+! include "COMMON.SPLITELE"
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscagrad_ele=0.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscagrad_ele=gamm*(6*gamm-6.0d0)/rlamb_ele
+ else
+ sscagrad_ele=0.0d0
+ endif
+ return
+ end function sscagrad_ele
+!!!!!!!!!! PBCSCALE
+ real(kind=8) function sscale2(r,r_cc,r_ll)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm,r_cc,r_ll
+ if(r.lt.r_cc-r_ll) then
+ sscale2=1.0d0
+ else if(r.le.r_cc.and.r.ge.r_cc-r_ll) then
+ gamm=(r-(r_cc-r_ll))/r_ll
+ sscale2=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale2=0d0
+ endif
+ return
+ end function sscale2
+
+ real(kind=8) function sscagrad2(r,r_cc,r_ll)
+ real(kind=8) :: r,gamm,r_cc,r_ll
+! include "COMMON.SPLITELE"
+ if(r.lt.r_cc-r_ll) then
+ sscagrad2=0.0d0
+ else if(r.le.r_cc.and.r.ge.r_cc-r_ll) then
+ gamm=(r-(r_cc-r_ll))/r_ll
+ sscagrad2=gamm*(6*gamm-6.0d0)/r_ll
+ else
+ sscagrad2=0.0d0
+ endif
+ return
+ end function sscagrad2
+
+ real(kind=8) function sscalelip(r)
+ real(kind=8) r,gamm
+ sscalelip=1.0d0+r*r*(2.0d0*r-3.0d0)
+ return
+ end function sscalelip
+!C-----------------------------------------------------------------------
+ real(kind=8) function sscagradlip(r)
+ real(kind=8) r,gamm
+ sscagradlip=r*(6.0d0*r-6.0d0)
+ return
+ end function sscagradlip
+
+!!!!!!!!!!!!!!!
+!-----------------------------------------------------------------------------
+ subroutine elj_long(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the LJ potential of interaction.
+!
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.TORSION'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CONTACTS'
+ real(kind=8),parameter :: accur=1.0d-10
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
+!el local variables
+ integer :: i,iint,j,k,itypi,itypi1,itypj
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,rij,sss,rrij,fac,eps0ij
+ real(kind=8) :: e1,e2,evdwij,evdw,sslipi,ssgradlipi,&
+ sslipj,ssgradlipj,aa,bb
+! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
+!d & 'iend=',iend(i,iint)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
+ if (sss.lt.1.0d0) then
+ rrij=1.0D0/rij
+ eps0ij=eps(itypi,itypj)
+ fac=rrij**expon2
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=e1+e2
+ evdw=evdw+(1.0d0-sss)*evdwij
+!
+! Calculate the components of the gradient in DC and X
+!
+ fac=-rrij*(e1+evdwij)*(1.0d0-sss)
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+!******************************************************************************
+!
+! N O T E !!!
+!
+! To save time, the factor of EXPON has been extracted from ALL components
+! of GVDWC and GRADX. Remember to multiply them by this factor before further
+! use!
+!
+!******************************************************************************
+ return
+ end subroutine elj_long
+!-----------------------------------------------------------------------------
+ subroutine elj_short(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the LJ potential of interaction.
+!
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.TORSION'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CONTACTS'
+ real(kind=8),parameter :: accur=1.0d-10
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
+!el local variables
+ integer :: i,iint,j,k,itypi,itypi1,itypj,num_conti
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,rij,sss,rrij,fac,eps0ij
+ real(kind=8) :: e1,e2,evdwij,evdw,sslipi,ssgradlipi,&
+ sslipj,ssgradlipj
+! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+! Change 12/1/95
+ num_conti=0
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
+!d & 'iend=',iend(i,iint)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+! Change 12/1/95 to calculate four-body interactions
+ rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
+ if (sss.gt.0.0d0) then
+ rrij=1.0D0/rij
+ eps0ij=eps(itypi,itypj)
+ fac=rrij**expon2
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=e1+e2
+ evdw=evdw+sss*evdwij
+!
+! Calculate the components of the gradient in DC and X
+!
+ fac=-rrij*(e1+evdwij)*sss
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+!******************************************************************************
+!
+! N O T E !!!
+!
+! To save time, the factor of EXPON has been extracted from ALL components
+! of GVDWC and GRADX. Remember to multiply them by this factor before further
+! use!
+!
+!******************************************************************************
+ return
+ end subroutine elj_short
+!-----------------------------------------------------------------------------
+ subroutine eljk_long(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the LJK potential of interaction.
+!
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
+ logical :: scheck
+!el local variables
+ integer :: i,iint,j,k,itypi,itypi1,itypj
+ real(kind=8) :: rrij,r_inv_ij,xj,yj,zj,xi,yi,zi,fac,evdw,&
+ fac_augm,e_augm,rij,sss,r_shift_inv,e1,e2,evdwij
+! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ r_inv_ij=dsqrt(rrij)
+ rij=1.0D0/r_inv_ij
+ sss=sscale(rij/sigma(itypi,itypj))
+ if (sss.lt.1.0d0) then
+ r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
+ fac=r_shift_inv**expon
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=e_augm+e1+e2
+!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
+!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
+!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
+!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
+ evdw=evdw+(1.0d0-sss)*evdwij
+!
+! Calculate the components of the gradient in DC and X
+!
+ fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ fac=fac*(1.0d0-sss)
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+ return
+ end subroutine eljk_long
+!-----------------------------------------------------------------------------
+ subroutine eljk_short(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the LJK potential of interaction.
+!
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
+ logical :: scheck
+!el local variables
+ integer :: i,iint,j,k,itypi,itypi1,itypj
+ real(kind=8) :: rrij,r_inv_ij,xj,yj,zj,xi,yi,zi,fac,evdw,&
+ fac_augm,e_augm,rij,sss,r_shift_inv,e1,e2,evdwij,&
+ sslipi,ssgradlipi,sslipj,ssgradlipj,aa,bb
+! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ r_inv_ij=dsqrt(rrij)
+ rij=1.0D0/r_inv_ij
+ sss=sscale(rij/sigma(itypi,itypj))
+ if (sss.gt.0.0d0) then
+ r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
+ fac=r_shift_inv**expon
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=e_augm+e1+e2
+!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
+!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
+!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
+!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
+ evdw=evdw+sss*evdwij
+!
+! Calculate the components of the gradient in DC and X
+!
+ fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ fac=fac*sss
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ enddo
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwc(j,i)=expon*gvdwc(j,i)
+ gvdwx(j,i)=expon*gvdwx(j,i)
+ enddo
+ enddo
+ return
+ end subroutine eljk_short
+!-----------------------------------------------------------------------------
+ subroutine ebp_long(evdw)
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the Berne-Pechukas potential of interaction.
+!
+ use calc_data
+! implicit real(kind=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'
+ use comm_srutu
+!el integer :: icall
+!el common /srutu/ icall
+! double precision rrsave(maxdim)
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj
+ real(kind=8) :: rrij,xi,yi,zi,fac,sslipi,ssgradlipi,&
+ sslipj,ssgradlipj,aa,bb
+ real(kind=8) :: sss,e1,e2,evdw,sigm,epsi
+ evdw=0.0D0
+! print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+! if (icall.eq.0) then
+! lprn=.true.
+! else
+ lprn=.false.
+! endif
+!el ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+!el ind=ind+1
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+!chi1=chi(itypi,itypj)
+!chi2=chi(itypj,itypi)
+!chi12=chi1*chi2
+!chip1=chip(itypi)
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.lt.1.0d0) then
+
+ ! Calculate the angle-dependent terms of energy & contributions to derivatives.
+ call sc_angular
+ ! Calculate whole angle-dependent part of epsilon and contributions
+ ! to its derivatives
+ fac=(rrij*sigsq)**expon2
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*(1.0d0-sss)
+ if (lprn) then
+ sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
+ !d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
+ !d & restyp(itypi,1),i,restyp(itypj,1),j,
+ !d & epsi,sigm,chi1,chi2,chip1,chip2,
+ !d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
+ !d & om1,om2,om12,1.0D0/dsqrt(rrij),
+ !d & evdwij
+ endif
+ ! Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)
+ sigder=fac/sigsq
+ fac=rrij*fac
+ ! Calculate radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ ! Calculate the angular part of the gradient and sum add the contributions
+ ! to the appropriate components of the Cartesian gradient.
+ call sc_grad_scale(1.0d0-sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ ! stop
+ return
+ end subroutine ebp_long
+ !-----------------------------------------------------------------------------
+ subroutine ebp_short(evdw)
+ !
+ ! This subroutine calculates the interaction energy of nonbonded side chains
+ ! assuming the Berne-Pechukas potential of interaction.
+ !
+ use calc_data
+! implicit real(kind=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'
+ use comm_srutu
+ !el integer :: icall
+ !el common /srutu/ icall
+! double precision rrsave(maxdim)
+ logical :: lprn
+ !el local variables
+ integer :: iint,itypi,itypi1,itypj
+ real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi
+ real(kind=8) :: sss,e1,e2,evdw,aa,bb, &
+ sslipi,ssgradlipi,sslipj,ssgradlipj
+ evdw=0.0D0
+ ! print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ ! if (icall.eq.0) then
+ ! lprn=.true.
+ ! else
+ lprn=.false.
+ ! endif
+ !el ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ ! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+ !
+ ! Calculate SC interaction energy.
+ !
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ !el ind=ind+1
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+ ! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.gt.0.0d0) then
+
+! Calculate the angle-dependent terms of energy & contributions to derivatives.
+ call sc_angular
+! Calculate whole angle-dependent part of epsilon and contributions
+! to its derivatives
+ fac=(rrij*sigsq)**expon2
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*sss
+ if (lprn) then
+ sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
+!d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
+!d & epsi,sigm,chi1,chi2,chip1,chip2,
+!d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
+!d & om1,om2,om12,1.0D0/dsqrt(rrij),
+!d & evdwij
+ endif
+! Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)
+ sigder=fac/sigsq
+ fac=rrij*fac
+! Calculate radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+! Calculate the angular part of the gradient and sum add the contributions
+! to the appropriate components of the Cartesian gradient.
+ call sc_grad_scale(sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+! stop
+ return
+ end subroutine ebp_short
+!-----------------------------------------------------------------------------
+ subroutine egb_long(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the Gay-Berne potential of interaction.
+!
+ use calc_data
+! implicit real(kind=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'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
+ real(kind=8) :: sss,e1,e2,evdw,sss_grad
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,fracinbuf,sslipi,sslipj,&
+ ssgradlipi,ssgradlipj
+
+
+ evdw=0.0D0
+!cccc energy_dec=.false.
+! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+! if (icall.eq.0) lprn=.false.
+!el ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+! write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
+! write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
+!
+! Calculate SC interaction energy.
+!
+ 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'
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,' ss'
+! do k=j+1,iend(i,iint)
+!C search over all next residues
+! if (dyn_ss_mask(k)) then
+!C check if they are cysteins
+!C write(iout,*) 'k=',k
+
+!c write(iout,*) "PRZED TRI", evdwij
+! evdwij_przed_tri=evdwij
+! call triple_ssbond_ene(i,j,k,evdwij)
+!c if(evdwij_przed_tri.ne.evdwij) then
+!c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+!c endif
+
+!c write(iout,*) "PO TRI", evdwij
+!C call the energy function that removes the artifical triple disulfide
+!C bond the soubroutine is located in ssMD.F
+! evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,'tss'
+! endif!dyn_ss_mask(k)
+! enddo! k
+
+ ELSE
+!el ind=ind+1
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+! & 1.0d0/vbld(j+nres)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
+ sig0ij=sigma(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+! Searching for nearest neighbour
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+ if (sss_ele_cut.le.0.0) cycle
+ if (sss.lt.1.0d0) then
+
+! Calculate angle-dependent terms of energy and contributions to their
+! derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+! for diagnostics; uncomment
+! rij_shift=1.2*sig0ij
+! I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
+!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
+ return
+ endif
+ sigder=-sig*sigsq
+!---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa
+ e2=fac*bb
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
+! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*(1.0d0-sss)*sss_ele_cut
+ if (lprn) then
+ sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,1),i,restyp(itypj,1),j,&
+ epsi,sigm,chi1,chi2,chip1,chip2,&
+ eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'evdw',i,j,evdwij
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,"egb_long"
+
+! Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ *rij-sss_grad/(1.0-sss)*rij &
+ /sigmaii(itypi,itypj))
+! fac=0.0d0
+! Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+! Calculate angular part of the gradient.
+ call sc_grad_scale(1.0d0-sss)
+ ENDIF !mask_dyn_ss
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+! write (iout,*) "Number of loop steps in EGB:",ind
+!ccc energy_dec=.false.
+ return
+ end subroutine egb_long
+!-----------------------------------------------------------------------------
+ subroutine egb_short(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the Gay-Berne potential of interaction.
+!
+ use calc_data
+! implicit real(kind=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'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap,countss
+ real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
+ real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,fracinbuf,sslipi,sslipj,&
+ ssgradlipi,ssgradlipj
+ evdw=0.0D0
+!cccc energy_dec=.false.
+! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+ countss=0
+! if (icall.eq.0) lprn=.false.
+!el ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+ 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
+ countss=countss+1
+ call dyn_ssbond_ene(i,j,evdwij,countss)
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,' ss'
+ do k=j+1,iend(i,iint)
+!C search over all next residues
+ if (dyn_ss_mask(k)) then
+!C check if they are cysteins
+!C write(iout,*) 'k=',k
+
+!c write(iout,*) "PRZED TRI", evdwij
+! evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
+!c if(evdwij_przed_tri.ne.evdwij) then
+!c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+!c endif
+
+!c write(iout,*) "PO TRI", evdwij
+!C call the energy function that removes the artifical triple disulfide
+!C bond the soubroutine is located in ssMD.F
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,'tss'
+ endif!dyn_ss_mask(k)
+ enddo! k
+ ELSE
+
+! typj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ dscj_inv=dsc_inv(itypj)
+! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+! & 1.0d0/vbld(j+nres)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
+ sig0ij=sigma(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+! Searching for nearest neighbour
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
+ if (sss_ele_cut.le.0.0) cycle
+
+ if (sss.gt.0.0d0) then
+
+! Calculate angle-dependent terms of energy and contributions to their
+! derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+! for diagnostics; uncomment
+! rij_shift=1.2*sig0ij
+! I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
+!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
+ return
+ endif
+ sigder=-sig*sigsq
+!---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa
+ e2=fac*bb
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
+! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*sss*sss_ele_cut
+ if (lprn) then
+ sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,1),i,restyp(itypj,1),j,&
+ epsi,sigm,chi1,chi2,chip1,chip2,&
+ eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'evdw',i,j,evdwij
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,"egb_short"
+
+! Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ *rij+sss_grad/sss*rij &
+ /sigmaii(itypi,itypj))
+
+! fac=0.0d0
+! Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+! Calculate angular part of the gradient.
+ call sc_grad_scale(sss)
+ endif
+ ENDIF !mask_dyn_ss
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+! write (iout,*) "Number of loop steps in EGB:",ind
+!ccc energy_dec=.false.
+ return
+ end subroutine egb_short
+!-----------------------------------------------------------------------------
+ subroutine egbv_long(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the Gay-Berne-Vorobjev potential of interaction.
+!
+ use calc_data
+! implicit real(kind=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'
+ use comm_srutu
+!el integer :: icall
+!el common /srutu/ icall
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj
+ real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,r0ij,sig,sig0ij,&
+ sslipi,ssgradlipi,sslipj,ssgradlipj,aa,bb
+ real(kind=8) :: sss,e1,e2,evdw,fac_augm,e_augm,rij_shift
+ evdw=0.0D0
+! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+! if (icall.eq.0) lprn=.true.
+!el ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+!el ind=ind+1
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma(itypi,itypj)
+ r0ij=r0(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.lt.1.0d0) then
+
+! Calculate angle-dependent terms of energy and contributions to their
+! derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+r0ij
+! I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+!---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
+ if (lprn) then
+ sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,1),i,restyp(itypj,1),j,&
+ epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
+ chi1,chi2,chip1,chip2,&
+ eps1,eps2rt**2,eps3rt**2,&
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij+e_augm
+ endif
+! Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac-2*expon*rrij*e_augm
+! Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+! Calculate angular part of the gradient.
+ call sc_grad_scale(1.0d0-sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ end subroutine egbv_long
+!-----------------------------------------------------------------------------
+ subroutine egbv_short(evdw)
+!
+! This subroutine calculates the interaction energy of nonbonded side chains
+! assuming the Gay-Berne-Vorobjev potential of interaction.
+!
+ use calc_data
+! implicit real(kind=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'
+ use comm_srutu
+!el integer :: icall
+!el common /srutu/ icall
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj
+ real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,rij_shift,&
+ sslipi,ssgradlipi, sslipj,ssgradlipj,aa,bb
+ real(kind=8) :: sss,e1,e2,evdw,r0ij,sig,sig0ij,fac_augm,e_augm
+ evdw=0.0D0
+! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+ evdw=0.0D0
+ lprn=.false.
+! if (icall.eq.0) lprn=.true.
+!el ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=itype(i,1)
+ if (itypi.eq.ntyp1) cycle
+ itypi1=itype(i+1,1)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+!el ind=ind+1
+ itypj=itype(j,1)
+ if (itypj.eq.ntyp1) cycle
+! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma(itypi,itypj)
+ r0ij=r0(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+
+ if (sss.gt.0.0d0) then
+
+! Calculate angle-dependent terms of energy and contributions to their
+! derivatives.
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+r0ij
+! I hate to put IF's in the loops, but here don't have another choice!!!!
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+!---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+(evdwij+e_augm)*sss
+ if (lprn) then
+ sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,1),i,restyp(itypj,1),j,&
+ epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
+ chi1,chi2,chip1,chip2,&
+ eps1,eps2rt**2,eps3rt**2,&
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij+e_augm
+ endif
+! Calculate gradient components.
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac-2*expon*rrij*e_augm
+! Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+! Calculate angular part of the gradient.
+ call sc_grad_scale(sss)
+ endif
+ enddo ! j
+ enddo ! iint
+ enddo ! i
+ end subroutine egbv_short
+!-----------------------------------------------------------------------------
+ subroutine eelec_scale(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+!
+! This subroutine calculates the average interaction energy and its gradient
+! in the virtual-bond vectors between non-adjacent peptide groups, based on
+! the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
+! The potential depends both on the distance of peptide-group centers and on
+! the orientation of the CA-CA virtual bonds.
+!
+! implicit real(kind=8) (a-h,o-z)
+
+ use comm_locel
+#ifdef MPI
+ include 'mpif.h'
+#endif
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SETUP'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+! include 'COMMON.TORSION'
+! include 'COMMON.VECTORS'
+! include 'COMMON.FFIELD'
+! include 'COMMON.TIME1'
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg
+ real(kind=8),dimension(3,3) ::erder,uryg,urzg,vryg,vrzg
+ real(kind=8),dimension(2,2) :: acipa !el,a_temp
+!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
+ real(kind=8),dimension(4) :: muij
+!el integer :: num_conti,j1,j2
+!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
+!el dz_normi,xmedi,ymedi,zmedi
+!el common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,&
+!el dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
+!el num_conti,j1,j2
+! 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+#ifdef MOMENT
+ real(kind=8) :: scal_el=1.0d0
+#else
+ real(kind=8) :: scal_el=0.5d0
+#endif
+! 12/13/98
+! 13-go grudnia roku pamietnego...
+ real(kind=8),dimension(3,3) :: unmat=reshape((/1.0d0,0.0d0,0.0d0,&
+ 0.0d0,1.0d0,0.0d0,&
+ 0.0d0,0.0d0,1.0d0/),shape(unmat))
+!el local variables
+ integer :: i,j,k
+ real(kind=8) :: fac
+ real(kind=8) :: dxj,dyj,dzj
+ real(kind=8) :: ees,evdw1,eel_loc,eello_turn3,eello_turn4
+
+! allocate(num_cont_hb(nres)) !(maxres)
+!d write(iout,*) 'In EELEC'
+!d do i=1,nloctyp
+!d write(iout,*) 'Type',i
+!d write(iout,*) 'B1',B1(:,i)
+!d write(iout,*) 'B2',B2(:,i)
+!d write(iout,*) 'CC',CC(:,:,i)
+!d write(iout,*) 'DD',DD(:,:,i)
+!d write(iout,*) 'EE',EE(:,:,i)
+!d enddo
+!d call check_vecgrad
+!d stop
+ if (icheckgrad.eq.1) then
+ do i=1,nres-1
+ fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
+ do k=1,3
+ dc_norm(k,i)=dc(k,i)*fac
+ enddo
+! write (iout,*) 'i',i,' fac',fac
+ enddo
+ endif
+ 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
+! call vec_and_deriv
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
+! print *, "before set matrices"
+ call set_matrices
+! print *,"after set catices"
+#ifdef TIMING
+ time_mat=time_mat+MPI_Wtime()-time01
+#endif
+ endif
+!d do i=1,nres-1
+!d write (iout,*) 'i=',i
+!d do k=1,3
+!d write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
+!d enddo
+!d do k=1,3
+!d write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
+!d & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
+!d enddo
+!d enddo
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ eel_loc=0.0d0
+ eello_turn3=0.0d0
+ eello_turn4=0.0d0
+!el ind=0
+ do i=1,nres
+ num_cont_hb(i)=0
+ enddo
+!d print '(a)','Enter EELEC'
+!d write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
+! if (.not.allocated(gel_loc_loc)) allocate(gel_loc_loc(nres)) !(maxvar)(maxvar=6*maxres)
+! if (.not.allocated(gcorr_loc)) allocate(gcorr_loc(nres)) !(maxvar)(maxvar=6*maxres)
+ do i=1,nres
+ gel_loc_loc(i)=0.0d0
+ gcorr_loc(i)=0.0d0
+ enddo
+!
+!
+! 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
+!
+! Loop over i,i+2 and i,i+3 pairs of the peptide groups
+!
+ do i=iturn3_start,iturn3_end
+ if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
+ num_conti=0
+ call eelecij_scale(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 (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+3,1).eq.ntyp1 &
+ .or. itype(i+4,1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+
+ call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
+
+ num_conti=num_cont_hb(i)
+ call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
+ if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
+ call eturn4(i,eello_turn4)
+ num_cont_hb(i)=num_conti
+ enddo ! i
+!
+! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
+!
+ do i=iatel_s,iatel_e
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
+! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+ num_conti=num_cont_hb(i)
+ do j=ielstart(i),ielend(i)
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
+ call eelecij_scale(i,j,ees,evdw1,eel_loc)
+ enddo ! j
+ num_cont_hb(i)=num_conti
+ enddo ! i
+! write (iout,*) "Number of loop steps in EELEC:",ind
+!d do i=1,nres
+!d write (iout,'(i3,3f10.5,5x,3f10.5)')
+!d & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
+!d enddo
+! 12/7/99 Adam eello_turn3 will be considered as a separate energy term
+!cc eel_loc=eel_loc+eello_turn3
+!d print *,"Processor",fg_rank," t_eelecij",t_eelecij
+ return
+ end subroutine eelec_scale
+!-----------------------------------------------------------------------------
+ subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
+! implicit real(kind=8) (a-h,o-z)
+
+ use comm_locel
+! include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+! include 'COMMON.CONTROL'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+! include 'COMMON.TORSION'
+! include 'COMMON.VECTORS'
+! include 'COMMON.FFIELD'
+! include 'COMMON.TIME1'
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg,xtemp
+ real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
+ real(kind=8),dimension(2,2) :: acipa !el,a_temp
+!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
+ real(kind=8),dimension(4) :: muij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,sss_grad
+ integer xshift,yshift,zshift
+
+!el integer :: num_conti,j1,j2
+!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
+!el dz_normi,xmedi,ymedi,zmedi
+!el common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,&
+!el dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
+!el num_conti,j1,j2
+! 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+#ifdef MOMENT
+ real(kind=8) :: scal_el=1.0d0
+#else
+ real(kind=8) :: scal_el=0.5d0
+#endif
+! 12/13/98
+! 13-go grudnia roku pamietnego...
+ real(kind=8),dimension(3,3) :: unmat=reshape((/1.0d0,0.0d0,0.0d0,&
+ 0.0d0,1.0d0,0.0d0,&
+ 0.0d0,0.0d0,1.0d0/),shape(unmat))
+!el local variables
+ integer :: i,j,k,l,iteli,itelj,kkk,kkll,m,isubchap
+ real(kind=8) :: aaa,bbb,ael6i,ael3i,dxj,dyj,dzj
+ real(kind=8) :: xj,yj,zj,rij,rrmij,rmij,sss,r3ij,r6ij,fac
+ real(kind=8) :: cosa,cosb,cosg,ev1,ev2,fac3,fac4,evdwij
+ real(kind=8) :: el1,el2,eesij,ees0ij,r0ij,fcont,fprimcont
+ real(kind=8) :: ees0tmp,ees0pij1,ees0mij1,ees0pijp,ees0mijp
+ real(kind=8) :: ees,evdw1,eel_loc,eel_loc_ij,dx_normj,dy_normj,&
+ dz_normj,facvdw,facel,fac1,facr,ecosa,ecosb,ecosg,&
+ ury,urz,vry,vrz,a22der,a23der,a32der,a33der,cosa4,&
+ wij,cosbg1,cosbg2,ees0pij,ees0mij,fac3p,ecosa1,ecosb1,&
+ ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,ecosgp,&
+ ecosam,ecosbm,ecosgm,ghalf,time00,faclipij,faclipij2
+! integer :: maxconts
+! maxconts = nres/4
+! allocate(gacontp_hb1(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(gacontp_hb2(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(gacontp_hb3(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(gacontm_hb1(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(gacontm_hb2(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(gacontm_hb3(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(gacont_hbr(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(grij_hb_cont(3,maxconts,nres)) !(3,maxconts,maxres) ! (maxconts=maxres/4)
+! allocate(facont_hb(maxconts,nres)) !(maxconts,maxres)
+! allocate(ees0p(maxconts,nres)) !(maxconts,maxres)
+! allocate(ees0m(maxconts,nres)) !(maxconts,maxres)
+! allocate(d_cont(maxconts,nres)) !(maxconts,maxres)
+! allocate(jcont_hb(maxconts,nres)) !(maxconts,maxres)
+
+! allocate(a_chuj(2,2,maxconts,nres)) !(2,2,maxconts,maxres)
+! allocate(a_chuj_der(2,2,3,5,maxconts,nres)) !(2,2,3,5,maxconts,maxres)
+
+#ifdef MPI
+ time00=MPI_Wtime()
+#endif
+!d write (iout,*) "eelecij",i,j
+!el ind=ind+1
+ iteli=itel(i)
+ itelj=itel(j)
+ if (j.eq.i+2 .and. itelj.eq.2) iteli=2
+ aaa=app(iteli,itelj)
+ bbb=bpp(iteli,itelj)
+ ael6i=ael6(iteli,itelj)
+ ael3i=ael3(iteli,itelj)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+! xj=c(1,j)+0.5D0*dxj-xmedi
+! yj=c(2,j)+0.5D0*dyj-ymedi
+! zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0
+ faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
+ rij=xj*xj+yj*yj+zj*zj
+ rrmij=1.0D0/rij
+ rij=dsqrt(rij)
+ rmij=1.0D0/rij
+! For extracting the short-range part of Evdwpp
+ sss=sscale(rij/rpp(iteli,itelj))
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+ sss_grad=sscale_grad((rij/rpp(iteli,itelj)))
+! sss_ele_cut=1.0d0
+! sss_ele_grad=0.0d0
+ if (sss_ele_cut.le.0.0) go to 128
+
+ r3ij=rrmij*rmij
+ r6ij=r3ij*r3ij
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ fac=cosa-3.0D0*cosb*cosg
+ ev1=aaa*r6ij*r6ij
+! 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
+ if (j.eq.i+2) ev1=scal_el*ev1
+ ev2=bbb*r6ij
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ evdwij=ev1+ev2
+ el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
+ el2=fac4*fac
+ eesij=el1+el2
+! 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*sss_ele_cut
+ evdw1=evdw1+evdwij*(1.0d0-sss)*sss_ele_cut
+!d write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
+!d & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
+!d & 1.0D0/dsqrt(rrmij),evdwij,eesij,
+!d & xmedi,ymedi,zmedi,xj,yj,zj
+
+ if (energy_dec) then
+ write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
+ write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+ endif
+
+!
+! Calculate contributions to the Cartesian gradient.
+!
+#ifdef SPLITELE
+ facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss_ele_cut
+ facel=-3*rrmij*(el1+eesij)*sss_ele_cut
+ fac1=fac
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
+!
+! Radial derivatives. First process both termini of the fragment (i,j)
+!
+ ggg(1)=facel*xj+sss_ele_grad*rmij*eesij*xj
+ ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj
+ ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj
+! do k=1,3
+! ghalf=0.5D0*ggg(k)
+! gelc(k,i)=gelc(k,i)+ghalf
+! gelc(k,j)=gelc(k,j)+ghalf
+! enddo
+! 9/28/08 AL Gradient compotents will be summed only at the end
+ do k=1,3
+ gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+ gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+ enddo
+!
+! Loop over residues i+1 thru j-1.
+!
+!grad do k=i+1,j-1
+!grad do l=1,3
+!grad gelc(l,k)=gelc(l,k)+ggg(l)
+!grad enddo
+!grad enddo
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj*(1.0d0-sss) &
+ -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*xj
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj*(1.0d0-sss) &
+ -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*yj
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj*(1.0d0-sss) &
+ -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*zj
+! do k=1,3
+! ghalf=0.5D0*ggg(k)
+! gvdwpp(k,i)=gvdwpp(k,i)+ghalf
+! gvdwpp(k,j)=gvdwpp(k,j)+ghalf
+! enddo
+! 9/28/08 AL Gradient compotents will be summed only at the end
+ do k=1,3
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ enddo
+!
+! Loop over residues i+1 thru j-1.
+!
+!grad do k=i+1,j-1
+!grad do l=1,3
+!grad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
+!grad enddo
+!grad enddo
+#else
+ facvdw=(ev1+evdwij)*(1.0d0-sss)*sss_ele_cut
+ facel=(el1+eesij)*sss_ele_cut
+ fac1=fac
+ fac=-3*rrmij*(facvdw+facvdw+facel)
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
+!
+! Radial derivatives. First process both termini of the fragment (i,j)
+!
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+! do k=1,3
+! ghalf=0.5D0*ggg(k)
+! gelc(k,i)=gelc(k,i)+ghalf
+! gelc(k,j)=gelc(k,j)+ghalf
+! enddo
+! 9/28/08 AL Gradient compotents will be summed only at the end
+ do k=1,3
+ gelc_long(k,j)=gelc(k,j)+ggg(k)
+ gelc_long(k,i)=gelc(k,i)-ggg(k)
+ enddo
+!
+! Loop over residues i+1 thru j-1.
+!
+!grad do k=i+1,j-1
+!grad do l=1,3
+!grad gelc(l,k)=gelc(l,k)+ggg(l)
+!grad enddo
+!grad enddo
+! 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
+ do k=1,3
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ enddo
+#endif
+!
+! Angular part
+!
+ ecosa=2.0D0*fac3*fac1+fac4
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
+ ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
+ do k=1,3
+ dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
+ dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
+ enddo
+!d print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
+!d & (dcosg(k),k=1,3)
+ do k=1,3
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k) )*sss_ele_cut
+ enddo
+! do k=1,3
+! ghalf=0.5D0*ggg(k)
+! gelc(k,i)=gelc(k,i)+ghalf
+! & +(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)+ghalf
+! & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+! & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+! enddo
+!grad do k=i+1,j-1
+!grad do l=1,3
+!grad gelc(l,k)=gelc(l,k)+ggg(l)
+!grad enddo
+!grad 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)&
+ *sss_ele_cut
+ gelc(k,j)=gelc(k,j) &
+ +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
+ + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+ *sss_ele_cut
+ gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+ gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+ enddo
+ 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
+!
+! 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
+! energy of a peptide unit is assumed in the form of a second-order
+! Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
+! Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
+! are computed for EVERY pair of non-contiguous peptide groups.
+!
+ if (j.lt.nres-1) then
+ j1=j+1
+ j2=j-1
+ else
+ j1=j-1
+ j2=j-2
+ endif
+ kkk=0
+ do k=1,2
+ do l=1,2
+ kkk=kkk+1
+ muij(kkk)=mu(k,i)*mu(l,j)
+ enddo
+ enddo
+!d write (iout,*) 'EELEC: i',i,' j',j
+!d write (iout,*) 'j',j,' j1',j1,' j2',j2
+!d write(iout,*) 'muij',muij
+ ury=scalar(uy(1,i),erij)
+ urz=scalar(uz(1,i),erij)
+ vry=scalar(uy(1,j),erij)
+ vrz=scalar(uz(1,j),erij)
+ a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
+ a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
+ a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
+ a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
+ fac=dsqrt(-ael6i)*r3ij
+ a22=a22*fac
+ a23=a23*fac
+ a32=a32*fac
+ a33=a33*fac
+!d write (iout,'(4i5,4f10.5)')
+!d & i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33
+!d write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
+!d write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
+!d & uy(:,j),uz(:,j)
+!d write (iout,'(4f10.5)')
+!d & scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
+!d & scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
+!d write (iout,'(4f10.5)') ury,urz,vry,vrz
+!d write (iout,'(9f10.5/)')
+!d & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
+! Derivatives of the elements of A in virtual-bond vectors
+ call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
+ do k=1,3
+ uryg(k,1)=scalar(erder(1,k),uy(1,i))
+ uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
+ uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
+ urzg(k,1)=scalar(erder(1,k),uz(1,i))
+ urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
+ urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
+ vryg(k,1)=scalar(erder(1,k),uy(1,j))
+ vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
+ vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
+ vrzg(k,1)=scalar(erder(1,k),uz(1,j))
+ vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
+ vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
+ enddo
+! Compute radial contributions to the gradient
+ facr=-3.0d0*rrmij
+ a22der=a22*facr
+ a23der=a23*facr
+ a32der=a32*facr
+ a33der=a33*facr
+ agg(1,1)=a22der*xj
+ agg(2,1)=a22der*yj
+ agg(3,1)=a22der*zj
+ agg(1,2)=a23der*xj
+ agg(2,2)=a23der*yj
+ agg(3,2)=a23der*zj
+ agg(1,3)=a32der*xj
+ agg(2,3)=a32der*yj
+ agg(3,3)=a32der*zj
+ agg(1,4)=a33der*xj
+ agg(2,4)=a33der*yj
+ agg(3,4)=a33der*zj
+! Add the contributions coming from er
+ fac3=-3.0d0*fac
+ do k=1,3
+ agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
+ agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
+ agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
+ agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
+ enddo
+ do k=1,3
+! Derivatives in DC(i)
+!grad ghalf1=0.5d0*agg(k,1)
+!grad ghalf2=0.5d0*agg(k,2)
+!grad ghalf3=0.5d0*agg(k,3)
+!grad ghalf4=0.5d0*agg(k,4)
+ aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j)) &
+ -3.0d0*uryg(k,2)*vry)!+ghalf1
+ aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j)) &
+ -3.0d0*uryg(k,2)*vrz)!+ghalf2
+ aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j)) &
+ -3.0d0*urzg(k,2)*vry)!+ghalf3
+ aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j)) &
+ -3.0d0*urzg(k,2)*vrz)!+ghalf4
+! Derivatives in DC(i+1)
+ aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j)) &
+ -3.0d0*uryg(k,3)*vry)!+agg(k,1)
+ aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j)) &
+ -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
+ aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j)) &
+ -3.0d0*urzg(k,3)*vry)!+agg(k,3)
+ aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j)) &
+ -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
+! Derivatives in DC(j)
+ aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i)) &
+ -3.0d0*vryg(k,2)*ury)!+ghalf1
+ aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i)) &
+ -3.0d0*vrzg(k,2)*ury)!+ghalf2
+ aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i)) &
+ -3.0d0*vryg(k,2)*urz)!+ghalf3
+ aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i)) &
+ -3.0d0*vrzg(k,2)*urz)!+ghalf4
+! Derivatives in DC(j+1) or DC(nres-1)
+ aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i)) &
+ -3.0d0*vryg(k,3)*ury)
+ aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i)) &
+ -3.0d0*vrzg(k,3)*ury)
+ aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i)) &
+ -3.0d0*vryg(k,3)*urz)
+ aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i)) &
+ -3.0d0*vrzg(k,3)*urz)
+!grad if (j.eq.nres-1 .and. i.lt.j-2) then
+!grad do l=1,4
+!grad aggj1(k,l)=aggj1(k,l)+agg(k,l)
+!grad enddo
+!grad endif
+ enddo
+ acipa(1,1)=a22
+ acipa(1,2)=a23
+ acipa(2,1)=a32
+ acipa(2,2)=a33
+ a22=-a22
+ a23=-a23
+ do l=1,2
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
+ enddo
+ enddo
+ if (j.lt.nres-1) then
+ a22=-a22
+ a32=-a32
+ do l=1,3,2
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
+ enddo
+ enddo
+ else
+ a22=-a22
+ a23=-a23
+ a32=-a32
+ a33=-a33
+ do l=1,4
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
+ enddo
+ enddo
+ endif
+ ENDIF ! WCORR
+ IF (wel_loc.gt.0.0d0) THEN
+! Contribution to the local-electrostatic energy coming from the i-j pair
+ eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) &
+ +a33*muij(4)
+! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+! print *,"EELLOC",i,gel_loc_loc(i-1)
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'eelloc',i,j,eel_loc_ij
+! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) !d
+
+ eel_loc=eel_loc+eel_loc_ij*sss_ele_cut
+! Partial derivatives in virtual-bond dihedral angles gamma
+ if (i.gt.1) &
+ gel_loc_loc(i-1)=gel_loc_loc(i-1)+ &
+ (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) &
+ +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) &
+ *sss_ele_cut
+ gel_loc_loc(j-1)=gel_loc_loc(j-1)+ &
+ (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) &
+ +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) &
+ *sss_ele_cut
+ xtemp(1)=xj
+ xtemp(2)=yj
+ xtemp(3)=zj
+
+! Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
+ do l=1,3
+ ggg(l)=(agg(l,1)*muij(1)+ &
+ agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))&
+ *sss_ele_cut &
+ +eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+
+ gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
+ gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
+!grad ghalf=0.5d0*ggg(l)
+!grad gel_loc(l,i)=gel_loc(l,i)+ghalf
+!grad gel_loc(l,j)=gel_loc(l,j)+ghalf
+ enddo
+!grad do k=i+1,j2
+!grad do l=1,3
+!grad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
+!grad enddo
+!grad enddo
+! 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))&
+ *sss_ele_cut
+
+ 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))&
+ *sss_ele_cut
+
+ 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))&
+ *sss_ele_cut
+
+ 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))&
+ *sss_ele_cut
+
+ enddo
+ ENDIF
+! Change 12/26/95 to calculate four-body contributions to H-bonding energy
+! if (j.gt.i+1 .and. num_conti.le.maxconts) then
+ if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0 &
+ .and. num_conti.le.maxconts) then
+! write (iout,*) i,j," entered corr"
+!
+! Calculate the contact function. The ith column of the array JCONT will
+! contain the numbers of atoms that make contacts with the atom I (of numbers
+! greater than I). The arrays FACONT and GACONT will contain the values of
+! the contact function and its derivative.
+! r0ij=1.02D0*rpp(iteli,itelj)
+! r0ij=1.11D0*rpp(iteli,itelj)
+ r0ij=2.20D0*rpp(iteli,itelj)
+! r0ij=1.55D0*rpp(iteli,itelj)
+ call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
+!elwrite(iout,*) "num_conti",num_conti, "maxconts",maxconts
+ if (fcont.gt.0.0D0) then
+ num_conti=num_conti+1
+ if (num_conti.gt.maxconts) then
+!elwrite(iout,*) "num_conti",num_conti, "maxconts",maxconts
+ write (iout,*) 'WARNING - max. # of contacts exceeded;',&
+ ' will skip next contacts for this conf.',num_conti
+ else
+ jcont_hb(num_conti,i)=j
+!d write (iout,*) "i",i," j",j," num_conti",num_conti,
+!d & " jcont_hb",jcont_hb(num_conti,i)
+ IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. &
+ wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
+! 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
+! terms.
+ d_cont(num_conti,i)=rij
+!d write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
+! --- Electrostatic-interaction matrix ---
+ a_chuj(1,1,num_conti,i)=a22
+ a_chuj(1,2,num_conti,i)=a23
+ a_chuj(2,1,num_conti,i)=a32
+ a_chuj(2,2,num_conti,i)=a33
+! --- Gradient of rij
+ do kkk=1,3
+ grij_hb_cont(kkk,num_conti,i)=erij(kkk)
+ enddo
+ kkll=0
+ do k=1,2
+ do l=1,2
+ kkll=kkll+1
+ do m=1,3
+ a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
+ a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
+ a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
+ a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
+ a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
+ enddo
+ enddo
+ enddo
+ ENDIF
+ IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN