+!---------------------------------------------------------------------------
+ subroutine ecat_nucl(ecation_nucl)
+ integer i,j,k,subchap,itmp,inum,itypi,itypj
+ real(kind=8) :: xi,yi,zi,xj,yj,zj
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_nucl,Evan1,Evan2,Ecav,Egb,wdip1,wdip2, &
+ wvan1,wvan2,wgbsig,wgbeps,wgbchi,wgbchip,wcav1,wcav2,wcav3,wcav4, &
+ wcavsig,wcavchi,v1m,v1dpdx,wh2o,wc,Edip,rcs2,invrcs6,invrcs8,invrcs12, &
+ invrcs14,rcb,rcb2,invrcb,invrcb2,invrcb4,invrcb6,cosinus,cos2,dcosdcatconst, &
+ dcosdcalpconst,dcosdcmconst,rcav,rcav11,rcav12,constcav1,constcav2, &
+ constgb1,constgb2,constdvan1,constdvan2,sgb,sgb6,sgb7,sgb12,sgb13, &
+ cavnum,cavdenom,invcavdenom2,dcavnumdcos,dcavnumdr,dcavdenomdcos, &
+ dcavdenomdr,sslipi,ssgradlipi,sslipj,ssgradlipj,aa,bb
+ real(kind=8),dimension(3) ::gg,r,dEtotalCm,dEtotalCalp,dEvan1Cm,&
+ dEvan2Cm,cm1,cm,vcat,vsug,v1,v2,dx,vcm,dEdipCm,dEdipCalp, &
+ dEvan1Calp,dEvan2Cat,dEvan2Calp,dEtotalCat,dEdipCat,dEvan1Cat,dcosdcat, &
+ dcosdcalp,dcosdcm,dEgbdCat,dEgbdCalp,dEgbdCm,dEcavdCat,dEcavdCalp, &
+ dEcavdCm
+ real(kind=8),dimension(14) :: vcatnuclprm
+ ecation_nucl=0.0d0
+ if (nres_molec(5).eq.0) return
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ if ((itype(i,2).eq.ntyp1_molec(2))) cycle ! leave dummy atoms
+ xi=(c(1,i+nres))
+ yi=(c(2,i+nres))
+ zi=(c(3,i+nres))
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ do k=1,3
+ cm1(k)=dc(k,i+nres)
+ enddo
+ do j=itmp+1,itmp+nres_molec(5)
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ 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)
+
+ dist_init=xj**2+yj**2+zj**2
+
+ itypi=itype(i,2)
+ itypj=itype(j,5)
+ do k=1,13
+ vcatnuclprm(k)=catnuclprm(k,itypi,itypj)
+ enddo
+ do k=1,3
+ vcm(k)=c(k,i+nres)
+ vsug(k)=c(k,i)
+ vcat(k)=c(k,j)
+ enddo
+ do k=1,3
+ dx(k) = vcat(k)-vcm(k)
+ enddo
+ do k=1,3
+ v1(k)=dc(k,i+nres)
+ v2(k)=(vcat(k)-vsug(k))
+ enddo
+ v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+ v1dpdx = v1(1)*dx(1)+v1(2)*dx(2)+v1(3)*dx(3)
+! The weights of the energy function calculated from
+!The quantum mechanical Gaussian simulations of potassium and sodium with deoxynucleosides
+ wh2o=78
+ wdip1 = vcatnuclprm(1)
+ wdip1 = wdip1/wh2o !w1
+ wdip2 = vcatnuclprm(2)
+ wdip2 = wdip2/wh2o !w2
+ wvan1 = vcatnuclprm(3)
+ wvan2 = vcatnuclprm(4) !pis1
+ wgbsig = vcatnuclprm(5) !sigma0
+ wgbeps = vcatnuclprm(6) !epsi0
+ wgbchi = vcatnuclprm(7) !chi1
+ wgbchip = vcatnuclprm(8) !chip1
+ wcavsig = vcatnuclprm(9) !sig
+ wcav1 = vcatnuclprm(10) !b1
+ wcav2 = vcatnuclprm(11) !b2
+ wcav3 = vcatnuclprm(12) !b3
+ wcav4 = vcatnuclprm(13) !b4
+ wcavchi = vcatnuclprm(14) !chis1
+ rcs2 = v2(1)**2+v2(2)**2+v2(3)**2
+ invrcs6 = 1/rcs2**3
+ invrcs8 = invrcs6/rcs2
+ invrcs12 = invrcs6**2
+ invrcs14 = invrcs12/rcs2
+ rcb2 = dx(1)**2+dx(2)**2+dx(3)**2
+ rcb = sqrt(rcb2)
+ invrcb = 1/rcb
+ invrcb2 = invrcb**2
+ invrcb4 = invrcb2**2
+ invrcb6 = invrcb4*invrcb2
+ cosinus = v1dpdx/(v1m*rcb)
+ cos2 = cosinus**2
+ dcosdcatconst = invrcb2/v1m
+ dcosdcalpconst = invrcb/v1m**2
+ dcosdcmconst = invrcb2/v1m**2
+ do k=1,3
+ dcosdcat(k) = (v1(k)*rcb-dx(k)*v1m*cosinus)*dcosdcatconst
+ dcosdcalp(k) = (v1(k)*rcb*cosinus-dx(k)*v1m)*dcosdcalpconst
+ dcosdcm(k) = ((dx(k)-v1(k))*v1m*rcb+ &
+ cosinus*(dx(k)*v1m**2-v1(k)*rcb2))*dcosdcmconst
+ enddo
+ rcav = rcb/wcavsig
+ rcav11 = rcav**11
+ rcav12 = rcav11*rcav
+ constcav1 = 1-wcavchi*cos2
+ constcav2 = sqrt(constcav1)
+ constgb1 = 1/sqrt(1-wgbchi*cos2)
+ constgb2 = wgbeps*(1-wgbchip*cos2)**2
+ constdvan1 = 12*wvan1*wvan2**12*invrcs14
+ constdvan2 = 6*wvan1*wvan2**6*invrcs8
+!----------------------------------------------------------------------------
+!Gay-Berne term
+!---------------------------------------------------------------------------
+ sgb = 1/(1-constgb1+(rcb/wgbsig))
+ sgb6 = sgb**6
+ sgb7 = sgb6*sgb
+ sgb12 = sgb6**2
+ sgb13 = sgb12*sgb
+ Egb = constgb2*(sgb12-sgb6)
+ do k=1,3
+ dEgbdCat(k) = -constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+ +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcat(k)
+ dEgbdCm(k) = constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+ +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcm(k)
+ dEgbdCalp(k) = (constgb1**3*constgb2*wgbchi*cosinus &
+ *(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcalp(k)
+ enddo
+!----------------------------------------------------------------------------
+!cavity term
+!---------------------------------------------------------------------------
+ cavnum = sqrt(rcav*constcav2)+wcav2*rcav*constcav2-wcav3
+ cavdenom = 1+wcav4*rcav12*constcav1**6
+ Ecav = wcav1*cavnum/cavdenom
+ invcavdenom2 = 1/cavdenom**2
+ dcavnumdcos = -wcavchi*cosinus/constcav2 &
+ *(sqrt(rcav/constcav2)/2+wcav2*rcav)
+ dcavnumdr = (0.5*sqrt(constcav2/rcav)+wcav2*constcav2)/wcavsig
+ dcavdenomdcos = -12*wcav4*wcavchi*rcav12*constcav1**5*cosinus
+ dcavdenomdr = 12*wcav4/wcavsig*rcav11*constcav1**6
+ do k=1,3
+ dEcavdCat(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcat(k)+(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+ dEcavdCm(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcm(k)-(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+ dEcavdCalp(k) = (dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcalp(k)*wcav1*invcavdenom2
+ enddo
+!----------------------------------------------------------------------------
+!van der Waals and dipole-charge interaction energy
+!---------------------------------------------------------------------------
+ Evan1 = wvan1*wvan2**12*invrcs12
+ do k=1,3
+ dEvan1Cat(k) = -v2(k)*constdvan1
+ dEvan1Cm(k) = 0.0d0
+ dEvan1Calp(k) = v2(k)*constdvan1
+ enddo
+ Evan2 = -wvan1*wvan2**6*invrcs6
+ do k=1,3
+ dEvan2Cat(k) = v2(k)*constdvan2
+ dEvan2Cm(k) = 0.0d0
+ dEvan2Calp(k) = -v2(k)*constdvan2
+ enddo
+ Edip = wdip1*cosinus*invrcb2-wdip2*(1-cos2)*invrcb4
+ do k=1,3
+ dEdipCat(k) = (-2*wdip1*cosinus*invrcb4 &
+ +4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+ +dcosdcat(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+ dEdipCm(k) = (2*wdip1*cosinus*invrcb4 &
+ -4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+ +dcosdcm(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+ dEdipCalp(k) = dcosdcalp(k)*(wdip1*invrcb2 &
+ +2*wdip2*cosinus*invrcb4)
+ enddo
+ if (energy_dec) write (iout,'(2i5,4(a6,f7.3))') i,j, &
+ ' E GB ',Egb,' ECav ',Ecav,' Evdw ',Evan1+Evan2,' Edip ',Edip
+ ecation_nucl=ecation_nucl+Ecav+Egb+Edip+Evan1+Evan2
+ do k=1,3
+ dEtotalCat(k) = dEcavdCat(k)+dEvan1Cat(k)+dEvan2Cat(k) &
+ +dEgbdCat(k)+dEdipCat(k)
+ dEtotalCm(k) = dEcavdCm(k)+dEvan1Cm(k)+dEvan2Cm(k) &
+ +dEgbdCm(k)+dEdipCm(k)
+ dEtotalCalp(k) = dEcavdCalp(k)+dEgbdCalp(k)+dEvan1Calp(k) &
+ +dEdipCalp(k)+dEvan2Calp(k)
+ enddo
+ do k=1,3
+ gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+ gradnuclcatx(k,i)=gradnuclcatx(k,i)+dEtotalCm(k)
+ gradnuclcat(k,i)=gradnuclcat(k,i)+gg(k)
+ gradnuclcat(k,j)=gradnuclcat(k,j)+dEtotalCat(k)
+ enddo
+ enddo !j
+ enddo !i
+ return
+ end subroutine ecat_nucl
+