end subroutine sc_grad
subroutine sc_grad_cat
-! implicit real*8 (a-h,o-z)
use calc_data
-! include 'DIMENSIONS'
-! include 'COMMON.CHAIN'
-! include 'COMMON.DERIV'
-! include 'COMMON.CALC'
-! include 'COMMON.IOUNITS'
real(kind=8), dimension(3) :: dcosom1,dcosom2
-! print *,"wchodze"
eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 &
+dCAVdOM1+ dGCLdOM1+ dPOLdOM1
eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 &
! eom2=0.0d0
! eom12=evdwij*eps1_om12
! end diagnostics
-! write (iout,*) "eps2der",eps2der," eps3der",eps3der,&
-! " sigder",sigder
-! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
-! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
-!C print *,sss_ele_cut,'in sc_grad'
do k=1,3
dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
!
! Calculate the components of the gradient in DC and X
!
-!grad do k=i,j-1
-!grad do l=1,3
-!grad gvdwc(l,k)=gvdwc(l,k)+gg(l)
-!grad enddo
-!grad enddo
do l=1,3
gradpepcat(l,i)=gradpepcat(l,i)-gg(l)
gradpepcat(l,j)=gradpepcat(l,j)+gg(l)
enddo
end subroutine sc_grad_cat
+ subroutine sc_grad_cat_pep
+ use calc_data
+ real(kind=8), dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 &
+ +dCAVdOM1+ dGCLdOM1+ dPOLdOM1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 &
+ +dCAVdOM2+ dGCLdOM2+ dPOLdOM2
+
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
+ -2.0D0*alf12*eps3der+sigder*sigsq_om12&
+ +dCAVdOM12+ dGCLdOM12
+! diagnostics only
+! eom1=0.0d0
+! eom2=0.0d0
+! eom12=evdwij*eps1_om12
+! end diagnostics
+
+ do k=1,3
+ dcosom1(k) = rij * (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwc_pepbase(k,i)= gvdwc_pepbase(k,i) +0.5*(- gg(k)) &
+ + (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+ *dsci_inv*2.0 &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_pepbase(k,i+1)= gvdwc_pepbase(k,i+1) +0.5*(- gg(k)) &
+ - (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i))) &
+ *dsci_inv*2.0 &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gradpepcat(k,j)=gradpepcat(k,j)+gg(k)
+ enddo
+ end subroutine sc_grad_cat_pep
#ifdef CRYST_THETA
!-----------------------------------------------------------------------------
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-4
+ aincr=1.0D-5
write(iout,*) 'Calling CHECK_ECARTINT.'
nf=0
icall=0
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-7
+ aincr=1.0D-6
write(iout,*) 'Calling CHECK_ECARTINT.',aincr
nf=0
icall=0
integer troll,jj,istate
real (kind=8) :: dcosom1(3),dcosom2(3)
- ecations_prot_amber=0.0D0
+ evdw=0.0D0
if (nres_molec(5).eq.0) return
eps_out=80.0d0
! sss_ele_cut=1.0d0
do i=1,4
itmp=itmp+nres_molec(i)
enddo
+ go to 17
! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
do i=ibond_start,ibond_end
! print *,"I am in EVDW",i
itypi=iabs(itype(i,1))
+
! if (i.ne.47) cycle
- if (itypi.eq.ntyp1) cycle
+ if ((itypi.eq.ntyp1).or.(itypi.eq.10)) cycle
itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
gg(k) = 0.0d0
ENDDO
!c! Compute head-head and head-tail energies for each state
- isel = iabs(Qi) + iabs(Qj)
+ isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj)
IF (isel.eq.0) THEN
!c! No charges - do nothing
eheadtail = 0.0d0
!c write (iout,*) "Number of loop steps in EGB:",ind
!c energy_dec=.false.
! print *,"EVDW KURW",evdw,nres
-
- return
- end subroutine ecats_prot_amber
-
-!---------------------------------------------------------------------------
-! old for Ca2+
- subroutine ecat_prot(ecation_prot)
-! use calc_data
-! use comm_momo
- integer i,j,k,subchap,itmp,inum
- real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
- r7,r4,ecationcation
- real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
- dist_init,dist_temp,ecation_prot,rcal,rocal, &
- Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
- catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
- wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
- costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
- Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
- rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
- opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
- opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
- Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
- ndiv,ndivi
- real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
- gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
- dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
- tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
- v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
- dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
- dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
- dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
- dEvan1Cat
- real(kind=8),dimension(6) :: vcatprm
- ecation_prot=0.0d0
-! first lets calculate interaction with peptide groups
- if (nres_molec(5).eq.0) return
- itmp=0
- do i=1,4
- itmp=itmp+nres_molec(i)
- enddo
-! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+!!! return
+ 17 continue
do i=ibond_start,ibond_end
-! cycle
- if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
- xi=0.5d0*(c(1,i)+c(1,i+1))
- yi=0.5d0*(c(2,i)+c(2,i+1))
- zi=0.5d0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
+
+! print *,"I am in EVDW",i
+ itypi=10 ! the peptide group parameters are for glicine
+
+! if (i.ne.47) cycle
+ if ((itype(i,1).eq.ntyp1).or.itype(i+1,1).eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1,1))
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ xi=dmod(xi,boxxsize)
if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
+ yi=dmod(yi,boxysize)
if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
+ zi=dmod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
-
+ dxi=dc_norm(1,i)
+ dyi=dc_norm(2,i)
+ dzi=dc_norm(3,i)
+ dsci_inv=vbld_inv(i+1)/2.0
do j=itmp+1,itmp+nres_molec(5)
-! print *,"WTF",itmp,j,i
-! all parameters were for Ca2+ to approximate single charge divide by two
- ndiv=1.0
- if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
- wconst=78*ndiv
- wdip =1.092777950857032D2
- wdip=wdip/wconst
- wmodquad=-2.174122713004870D4
- wmodquad=wmodquad/wconst
- wquad1 = 3.901232068562804D1
- wquad1=wquad1/wconst
- wquad2 = 3
- wquad2=wquad2/wconst
- wvan1 = 0.1
- wvan2 = 6
-! itmp=0
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,5))
+ if ((itypj.eq.ntyp1)) cycle
+ CALL elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+ dscj_inv=0.0
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- xj=dmod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=dmod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=dmod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
xj=xj_safe+xshift*boxxsize
yj=yj_safe+yshift*boxysize
zj=zj_safe+zshift*boxzsize
zj_temp=zj
subchap=1
endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
xj=xj_temp-xi
yj=yj_temp-yi
zj=zj_temp-zi
- else
+ else
xj=xj_safe-xi
yj=yj_safe-yi
zj=zj_safe-zi
- endif
-! enddo
-! enddo
- rcpm = sqrt(xj**2+yj**2+zj**2)
- drcp_norm(1)=xj/rcpm
- drcp_norm(2)=yj/rcpm
- drcp_norm(3)=zj/rcpm
- dcmag=0.0
- do k=1,3
- dcmag=dcmag+dc(k,i)**2
- enddo
- dcmag=dsqrt(dcmag)
- do k=1,3
- myd_norm(k)=dc(k,i)/dcmag
- enddo
- costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+&
- drcp_norm(3)*myd_norm(3)
- rsecp = rcpm**2
- Ir = 1.0d0/rcpm
- Irsecp = 1.0d0/rsecp
- Irthrp = Irsecp/rcpm
- Irfourp = Irthrp/rcpm
- Irfiftp = Irfourp/rcpm
- Irsistp=Irfiftp/rcpm
- Irseven=Irsistp/rcpm
- Irtwelv=Irsistp*Irsistp
- Irthir=Irtwelv/rcpm
- sin2thet = (1-costhet*costhet)
- sinthet=sqrt(sin2thet)
- E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)&
- *sin2thet
- E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
- 2*wvan2**6*Irsistp)
- ecation_prot = ecation_prot+E1+E2
-! print *,"ecatprot",i,j,ecation_prot,rcpm
- dE1dr = -2*costhet*wdip*Irthrp-&
- (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
- dE2dr = 3*wquad1*wquad2*Irfourp- &
- 12*wvan1*wvan2**6*(wvan2**6*Irthir-Irseven)
- dEdcos = wdip*Irsecp-2*(wmodquad*Irfourp+wquad1*Irthrp)*costhet
- do k=1,3
- drdpep(k) = -drcp_norm(k)
- dcosdpep(k) = Ir*(costhet*drcp_norm(k)-myd_norm(k))
- dcosddci(k) = drcp_norm(k)/dcmag-costhet*myd_norm(k)/dcmag
- dEdpep(k) = (dE1dr+dE2dr)*drdpep(k)+dEdcos*dcosdpep(k)
- dEddci(k) = dEdcos*dcosddci(k)
- enddo
- do k=1,3
- gradpepcat(k,i)=gradpepcat(k,i)+0.5D0*dEdpep(k)-dEddci(k)
- gradpepcat(k,i+1)=gradpepcat(k,i+1)+0.5D0*dEdpep(k)+dEddci(k)
- gradpepcat(k,j)=gradpepcat(k,j)-dEdpep(k)
- enddo
- enddo ! j
- enddo ! i
-!------------------------------------------sidechains
-! do i=1,nres_molec(1)
- do i=ibond_start,ibond_end
- if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms
-! cycle
-! print *,i,ecation_prot
- xi=(c(1,i+nres))
- yi=(c(2,i+nres))
- zi=(c(3,i+nres))
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
- do k=1,3
- cm1(k)=dc(k,i+nres)
- enddo
- cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
- do j=itmp+1,itmp+nres_molec(5)
- ndiv=1.0
- if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+ endif
- xj=c(1,j)
- yj=c(2,j)
- zj=c(3,j)
- xj=dmod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=dmod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
+ dxj = 0.0d0! dc_norm( 1, nres+j )
+ dyj = 0.0d0!dc_norm( 2, nres+j )
+ dzj = 0.0d0! dc_norm( 3, nres+j )
+
+ itypi = 10
+ itypj = itype(j,5)
+! Parameters from fitting the analitical expressions to the PMF obtained by umbrella
+! sampling performed with amber package
+! alf1 = 0.0d0
+! alf2 = 0.0d0
+! alf12 = 0.0d0
+! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+ chi1 = chicat(itypi,itypj)
+ chis1 = chiscat(itypi,itypj)
+ chip1 = chippcat(itypi,itypj)
+! chi1=0.0d0
+! chis1=0.0d0
+! chip1=0.0d0
+ chi2=0.0
+ chip2=0.0
+ chis2=0.0
+! chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1cat(itypi,itypj)
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabcat(itypi,itypj)
+ if (eps_in.eq.0.0) eps_in=1.0
+
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! Rtail = 0.0d0
+
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0
+ ctail(k,2)=c(k,j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail location and distance calculations
+! dhead1
+ d1 = dheadcat(1, 1, itypi, itypj)
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
+ chead(k,2) = c(k, j)
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = Rtail - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq_cat(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_cat(itypi,itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij
+!#endif
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ print *,"TUT2",fac,chis1,sqom1,pom
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+ Chif = Rtail * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1cav * ( eagle + b2cav * ChiLambf - b3cav )
+ bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+
+ dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
+ dbot = 12.0d0 * b4cav * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+
+ dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
+ dbot = 12.0d0 * b4cav * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ DO k= 1, 3
+ ertail(k) = Rtail_distance(k)/Rtail
+ END DO
+ erdxi = scalar( ertail(1), dC_norm(1,i) )
+ erdxj = scalar( ertail(1), dC_norm(1,j) )
+ facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i)
+ facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - (( dFdR + gg(k) ) * pom)
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+ gg(k) = 0.0d0
+ ENDDO
+!c! Compute head-head and head-tail energies for each state
+ isel = 3
+!c! Dipole-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+ CALL edq_cat_pep(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! print *,"i,",i,eheadtail
+ eheadtail = 0.0d0
+
+ evdw = evdw + Fcav + eheadtail
+
+ IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+ Equad,evdwij+Fcav+eheadtail,evdw
+! evdw = evdw + Fcav + eheadtail
+
+! iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad_cat_pep
+! END IF
+!c!-------------------------------------------------------------------
+!c! NAPISY KONCOWE
+ END DO ! j
+ END DO ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+
+
+ return
+ end subroutine ecats_prot_amber
+
+!---------------------------------------------------------------------------
+! old for Ca2+
+ subroutine ecat_prot(ecation_prot)
+! use calc_data
+! use comm_momo
+ integer i,j,k,subchap,itmp,inum
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
+ r7,r4,ecationcation
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_prot,rcal,rocal, &
+ Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
+ catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
+ wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
+ costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
+ Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
+ rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
+ opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
+ opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
+ Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+ ndiv,ndivi
+ real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
+ gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
+ dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
+ tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
+ v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
+ dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
+ dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
+ dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
+ dEvan1Cat
+ real(kind=8),dimension(6) :: vcatprm
+ ecation_prot=0.0d0
+! first lets calculate interaction with peptide groups
+ if (nres_molec(5).eq.0) return
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+ do i=ibond_start,ibond_end
+! cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
+ xi=0.5d0*(c(1,i)+c(1,i+1))
+ yi=0.5d0*(c(2,i)+c(2,i+1))
+ zi=0.5d0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do j=itmp+1,itmp+nres_molec(5)
+! print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+ wconst=78*ndiv
+ wdip =1.092777950857032D2
+ wdip=wdip/wconst
+ wmodquad=-2.174122713004870D4
+ wmodquad=wmodquad/wconst
+ wquad1 = 3.901232068562804D1
+ wquad1=wquad1/wconst
+ wquad2 = 3
+ wquad2=wquad2/wconst
+ wvan1 = 0.1
+ wvan2 = 6
+! itmp=0
+
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+! enddo
+! enddo
+ rcpm = sqrt(xj**2+yj**2+zj**2)
+ drcp_norm(1)=xj/rcpm
+ drcp_norm(2)=yj/rcpm
+ drcp_norm(3)=zj/rcpm
+ dcmag=0.0
+ do k=1,3
+ dcmag=dcmag+dc(k,i)**2
+ enddo
+ dcmag=dsqrt(dcmag)
+ do k=1,3
+ myd_norm(k)=dc(k,i)/dcmag
+ enddo
+ costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+&
+ drcp_norm(3)*myd_norm(3)
+ rsecp = rcpm**2
+ Ir = 1.0d0/rcpm
+ Irsecp = 1.0d0/rsecp
+ Irthrp = Irsecp/rcpm
+ Irfourp = Irthrp/rcpm
+ Irfiftp = Irfourp/rcpm
+ Irsistp=Irfiftp/rcpm
+ Irseven=Irsistp/rcpm
+ Irtwelv=Irsistp*Irsistp
+ Irthir=Irtwelv/rcpm
+ sin2thet = (1-costhet*costhet)
+ sinthet=sqrt(sin2thet)
+ E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)&
+ *sin2thet
+ E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
+ 2*wvan2**6*Irsistp)
+ ecation_prot = ecation_prot+E1+E2
+! print *,"ecatprot",i,j,ecation_prot,rcpm
+ dE1dr = -2*costhet*wdip*Irthrp-&
+ (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
+ dE2dr = 3*wquad1*wquad2*Irfourp- &
+ 12*wvan1*wvan2**6*(wvan2**6*Irthir-Irseven)
+ dEdcos = wdip*Irsecp-2*(wmodquad*Irfourp+wquad1*Irthrp)*costhet
+ do k=1,3
+ drdpep(k) = -drcp_norm(k)
+ dcosdpep(k) = Ir*(costhet*drcp_norm(k)-myd_norm(k))
+ dcosddci(k) = drcp_norm(k)/dcmag-costhet*myd_norm(k)/dcmag
+ dEdpep(k) = (dE1dr+dE2dr)*drdpep(k)+dEdcos*dcosdpep(k)
+ dEddci(k) = dEdcos*dcosddci(k)
+ enddo
+ do k=1,3
+ gradpepcat(k,i)=gradpepcat(k,i)+0.5D0*dEdpep(k)-dEddci(k)
+ gradpepcat(k,i+1)=gradpepcat(k,i+1)+0.5D0*dEdpep(k)+dEddci(k)
+ gradpepcat(k,j)=gradpepcat(k,j)-dEdpep(k)
+ enddo
+ enddo ! j
+ enddo ! i
+!------------------------------------------sidechains
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms
+! cycle
+! print *,i,ecation_prot
+ xi=(c(1,i+nres))
+ yi=(c(2,i+nres))
+ zi=(c(3,i+nres))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do k=1,3
+ cm1(k)=dc(k,i+nres)
+ enddo
+ cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
+ do j=itmp+1,itmp+nres_molec(5)
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
zj=dmod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- dGLJdR * pom
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
- gradpepcatx(k,j) = gradpepcatx(k,j)+ dGCLdR * pom&
- + dGGBdR * pom+ dGCVdR * pom&
- + dPOLdR1 * (erhead_tail(k,1)&
- -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
- + dPOLdR2 * condor + dGLJdR * pom
+! gradpepcatx(k,j) = gradpepcatx(k,j)+ dGCLdR * pom&
+! + dGGBdR * pom+ dGCVdR * pom&
+! + dPOLdR1 * (erhead_tail(k,1)&
+! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
+! + dPOLdR2 * condor + dGLJdR * pom
gradpepcat(k,i) = gradpepcat(k,i) &
- dGCLdR * erhead(k)&
gradpepcatx(k,i) = gradpepcatx(k,i) &
- dPOLdR2 * (erhead_tail(k,2) &
-facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
- gradpepcatx(k,j) = gradpepcatx(k,j) &
- + dPOLdR2 * condor
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dPOLdR2 * condor
gradpepcat(k,i) = gradpepcat(k,i) &
- dPOLdR2 * erhead_tail(k,2)
!c!-------------------------------------------------------------------
!c! ecl
- sparrow = w1 * Qi * om1
- hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
ECL = sparrow / Rhead**2.0d0 &
- hawk / Rhead**4.0d0
!c!-------------------------------------------------------------------
dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ 4.0d0 * hawk / Rhead**5.0d0
!c! dF/dom1
- dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
!c! dF/dom2
- dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
!c--------------------------------------------------------------------
!c Polarization energy
!c Epol
!c!-------------------------------------------------------------------
!c! ecl
- sparrow = w1 * Qi * om1
- hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
ECL = sparrow / Rhead**2.0d0 &
- hawk / Rhead**4.0d0
!c!-------------------------------------------------------------------
dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ 4.0d0 * hawk / Rhead**5.0d0
!c! dF/dom1
- dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
!c! dF/dom2
- dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
!c--------------------------------------------------------------------
!c--------------------------------------------------------------------
!c Polarization energy
- dGLJdR * pom
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
- gradpepcatx(k,j) = gradpepcatx(k,j) &
- + dGCLdR * pom &
- + dPOLdR2 * condor &
- + dGLJdR * pom
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
gradpepcat(k,i) = gradpepcat(k,i) &
RETURN
END SUBROUTINE edq_cat
+ SUBROUTINE edq_cat_pep(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolcat(itypj,itypi)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+ print *,"CO?!.", w1,w2,Qj,om1
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c!-------------------------------------------------------------------
+
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+1)/2.0
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - dGCLdR * pom &
+! - dPOLdR2 * (erhead_tail(k,2) &
+! -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+! - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+
+ gradpepcat(k,i) = gradpepcat(k,i) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))
+
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_cat_pep
SUBROUTINE edd(ECL)
! IMPLICIT NONE
!c! Gay-berne var's
sig0ij = sigmacat( itypi,itypj )
chi1 = chicat( itypi, itypj )
-! chi2 = chi( itypj, itypi )
chi2 = 0.0d0
-! chi12 = chi1 * chi2
chi12 = 0.0d0
chip1 = chippcat( itypi, itypj )
-! chip2 = chipp( itypj, itypi )
chip2 = 0.0d0
-! chip12 = chip1 * chip2
chip12 = 0.0d0
-! chi1=0.0
-! chi2=0.0
-! chi12=0.0
-! chip1=0.0
-! chip2=0.0
-! chip12=0.0
!c! not used by momo potential, but needed by sc_angular which is shared
!c! by all energy_potential subroutines
alf1 = 0.0d0
alf2 = 0.0d0
alf12 = 0.0d0
-!c! location, location, location
-! xj = c( 1, nres+j ) - xi
-! yj = c( 2, nres+j ) - yi
-! zj = c( 3, nres+j ) - zi
dxj = dc_norm( 1, nres+j )
dyj = dc_norm( 2, nres+j )
dzj = dc_norm( 3, nres+j )
Qij = Qi * Qj
!c! chis1,2,12
chis1 = chiscat(itypi,itypj)
-! chis2 = chis(itypj,itypi)
chis2 = 0.0d0
-! chis12 = chis1 * chis2
chis12 = 0.0d0
sig1 = sigmap1cat(itypi,itypj)
sig2 = sigmap2cat(itypi,itypj)
!c! alpha factors from Fcav/Gcav
b1cav = alphasurcat(1,itypi,itypj)
-! b1cav=0.0
b2cav = alphasurcat(2,itypi,itypj)
b3cav = alphasurcat(3,itypi,itypj)
b4cav = alphasurcat(4,itypi,itypj)
Rtail = 0.0d0
DO k = 1, 3
ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i)
- ctail(k,2)=c(k,j+nres)-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
END DO
!c! tail distances will be themselves usefull elswhere
!c1 (in Gcav, for example)
!c! and moving by a d1 * dc_norm vector
!c! see unres publications for very informative images
chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
- chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+ chead(k,2) = c(k, j)
!c! distance
!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
RETURN
END SUBROUTINE elgrad_init_cat
+ SUBROUTINE elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+ use comm_momo
+ use calc_data
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = 10
+ itypj = itype(j,5)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigmacat( itypi,itypj )
+ chi1 = chicat( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chippcat( itypi, itypj )
+ chip2 = 0.0d0
+ chip12 = 0.0d0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ dxj = 0.0d0 !dc_norm( 1, nres+j )
+ dyj = 0.0d0 !dc_norm( 2, nres+j )
+ dzj = 0.0d0 !dc_norm( 3, nres+j )
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborncat(itypi,itypj) * rborncat(itypj,itypi)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = 0
+ Qj = ichargecat(itypj)
+! Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chiscat(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1cat(itypi,itypj)
+ sig2 = sigmap2cat(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+ wqd = wquadcat(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabcat(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0-dtailcat(1,itypi,itypj)*dc_norm(k,i)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publications for very informative images
+ chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
+ chead(k,2) = c(k, j)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!c!-------------------------------------------------------------------
+!c! zero everything that should be zero'ed
+ Egb = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ RETURN
+ END SUBROUTINE elgrad_init_cat_pep
double precision function tschebyshev(m,n,x,y)
implicit none
projects / unres4.git / commitdiff