! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
!----------------------------------------
real(kind=8),dimension(:,:),allocatable ::gradlipelec,gradlipbond,&
- gradlipang,gradliplj
+ gradlipang,gradliplj,gradpepmart, gradpepmartx
real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
! energies for protein nucleic acid interaction
real(kind=8) :: escbase,epepbase,escpho,epeppho
! energies for MARTINI
- real(kind=8) :: elipbond,elipang,elipelec,eliplj
+ real(kind=8) :: elipbond,elipang,elipelec,eliplj,elipidprot
#ifdef MPI
real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
weights_(49)=wpeppho
weights_(50)=wcatnucl
weights_(56)=wcat_tran
-
+ weights_(58)=wlip_prot
+ weights_(52)=wmartini
! wcatcat= weights(41)
! wcatprot=weights(42)
wscpho=weights(48)
wpeppho=weights(49)
wcatnucl=weights(50)
+ wmartini=weights(52)
wcat_tran=weights(56)
-
+ wlip_prot=weights(58)
! welpsb=weights(28)*fact(1)
!
! wcorr_nucl= weights(37)*fact(1)
! write (iout,*) "after make_SCp_inter_list"
if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
! write (iout,*) "after make_SCSC_inter_list"
-
+ if (nres_molec(4).gt.0) then
+ if (mod(itime_mat,imatupdate).eq.0) call make_lip_pep_list
+ endif
if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
if (nres_molec(5).gt.0) then
if (mod(itime_mat,imatupdate).eq.0) then
endif
call lipid_LJ(eliplj)
call lipid_elec(elipelec)
+ if (nres_molec(1).gt.0) then
+ call elip_prot(elipidprot)
+ else
+ elipidprot=0.0d0
+ endif
else
elipbond=0.0d0
elipang=0.0d0
energia(55)=elipelec
energia(56)=ecat_prottran
energia(57)=ecation_protang
+ energia(58)=elipidprot
! write(iout,*) elipelec,"elipelec"
! write(iout,*) elipang,"elipang"
! write(iout,*) eliplj,"eliplj"
ecation_nucl,ecat_prottran,ecation_protang
real(kind=8) :: escbase,epepbase,escpho,epeppho
integer :: i
- real(kind=8) :: elipbond,elipang,eliplj,elipelec
+ real(kind=8) :: elipbond,elipang,eliplj,elipelec,elipidprot
#ifdef MPI
integer :: ierr
real(kind=8) :: time00
elipelec=energia(55)
ecat_prottran=energia(56)
ecation_protang=energia(57)
+ elipidprot=energia(58)
! ecations_prot_amber=energia(50)
! energia(41)=ecation_prot
+wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
+wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
- +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +(elipbond+elipang+eliplj+elipelec)*wmartini&
+ +wcat_tran*ecat_prottran+ecation_protang&
+ +wlip_prot*elipidprot&
#ifdef WHAM_RUN
+0.0d0
#else
+wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
+wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
- +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +(elipbond+elipang+eliplj+elipelec)*wmartini&
+ +wcat_tran*ecat_prottran+ecation_protang&
+ +wlip_prot*elipidprot&
#ifdef WHAM_RUN
+0.0d0
#else
real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,&
ecation_nucl,ecat_prottran,ecation_protang
real(kind=8) :: escbase,epepbase,escpho,epeppho
- real(kind=8) :: elipbond,elipang,eliplj,elipelec
+ real(kind=8) :: elipbond,elipang,eliplj,elipelec,elipidprot
etot=energia(0)
evdw=energia(1)
evdw2=energia(2)
ecat_prottran=energia(56)
ecation_protang=energia(57)
ehomology_constr=energia(51)
-
+ elipidprot=energia(58)
! ecations_prot_amber=energia(50)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
ecationcation,wcatcat, &
escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
ecation_nucl,wcatnucl,ehomology_constr,&
- elipbond,elipang,eliplj,elipelec,etot
+ elipbond,elipang,eliplj,elipelec,elipidprot,wlip_prot,etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
'ELIPANG=',1pE16.6,'(matrini angle energy)'/&
'ELIPLJ=',1pE16.6,'(matrini Lennard-Jones energy)'/&
'ELIPELEC=',1pE16.6,'(matrini electrostatic energy)'/&
+ 'ELIPPROT=',1pE16.6,' WEIGHT=',1pD16.6,'(lipid prot)'/ &
'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, &
escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
- ecation_nucl,wcatnucl,ehomology_constr,etot
+ ecation_nucl,wcatnucl,ehomology_constr,elipidprot,wlip_prot,etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
'ELIPANG=',1pE16.6,'(matrini angle energy)'/&
'ELIPLJ=',1pE16.6,'(matrini Lennard-Jones energy)'/&
'ELIPELEC=',1pE16.6,'(matrini electrostatic energy)'/&
+ 'ELIPPROT=',1pE16.6,' WEIGHT=',1pD16.6,'(lipid prot)'/ &
'ETOT= ',1pE16.6,' (total)')
#endif
return
wpepbase*gvdwc_pepbase(j,i)+&
wscpho*gvdwc_scpho(j,i)+ &
wpeppho*gvdwc_peppho(j,i)+wcatnucl*gradnuclcat(j,i)+ &
- gradlipbond(j,i)+gradlipang(j,i)+gradliplj(j,i)+gradlipelec(j,i)+&
- wcat_tran*gradcattranc(j,i)+gradcatangc(j,i)
+ wmartini*(gradlipbond(j,i)+gradlipang(j,i)+gradliplj(j,i)+gradlipelec(j,i))+&
+ wcat_tran*gradcattranc(j,i)+gradcatangc(j,i)+&
+ wlip_prot*gradpepmart(j,i)
+
wpepbase*gvdwc_pepbase(j,i)+&
wscpho*gvdwc_scpho(j,i)+&
wpeppho*gvdwc_peppho(j,i)+wcatnucl*gradnuclcat(j,i)+&
- gradlipbond(j,i)+gradlipang(j,i)+gradliplj(j,i)+gradlipelec(j,i)+&
- wcat_tran*gradcattranc(j,i)+gradcatangc(j,i)
+ wmartini*(gradlipbond(j,i)+gradlipang(j,i)+gradliplj(j,i)+gradlipelec(j,i))+&
+ wcat_tran*gradcattranc(j,i)+gradcatangc(j,i)+&
+ wlip_prot*gradpepmart(j,i)
+wscbase*gvdwx_scbase(j,i) &
+wpepbase*gvdwx_pepbase(j,i)&
+wscpho*gvdwx_scpho(j,i)+wcatnucl*gradnuclcatx(j,i)&
- +wcat_tran*gradcattranx(j,i)+gradcatangx(j,i)
+ +wcat_tran*gradcattranx(j,i)+gradcatangx(j,i)&
+ +wlip_prot*gradpepmartx(j,i)
+
! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i)
enddo
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-5
+ aincr=graddelta
write(iout,*) 'Calling CHECK_ECARTINT.,kupa'
nf=0
icall=0
#endif
! print *, "before set matrices"
call set_matrices
-! print *,"after set martices"
+! print *,"after set catices"
#ifdef TIMING
time_mat=time_mat+MPI_Wtime()-time01
#endif
gradcattranx(j,i)=0.0d0
gradcatangx(j,i)=0.0d0
gradcatangc(j,i)=0.0d0
+ gradpepmart(j,i)=0.0d0
+ gradpepmartx(j,i)=0.0d0
duscdiff(j,i)=0.0d0
duscdiffx(j,i)=0.0d0
enddo
allocate(gvdwpp_nucl(3,-1:nres))
allocate(gradpepcat(3,-1:nres))
allocate(gradpepcatx(3,-1:nres))
+ allocate(gradpepmart(3,-1:nres))
+ allocate(gradpepmartx(3,-1:nres))
allocate(gradcatcat(3,-1:nres))
allocate(gradnuclcat(3,-1:nres))
allocate(gradnuclcatx(3,-1:nres))
allocate(newcontlistppj(300*nres))
allocate(newcontlistscpj(350*nres))
allocate(newcontlistj(300*nres))
+ allocate(newcontlistmartpi(300*nres))
+ allocate(newcontlistmartpj(300*nres))
+ allocate(newcontlistmartsci(300*nres))
+ allocate(newcontlistmartscj(300*nres))
+
allocate(newcontlistcatsctrani(300*nres))
allocate(newcontlistcatsctranj(300*nres))
allocate(newcontlistcatptrani(300*nres))
allocate(newcontlistcatpnormj(300*nres))
allocate(newcontlistcatcatnormi(900*nres))
allocate(newcontlistcatcatnormj(900*nres))
-
+
allocate(newcontlistcatscangi(300*nres))
allocate(newcontlistcatscangj(300*nres))
allocate(newcontlistcatscangfi(300*nres))
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
+ Rreal(1)=xj
+ Rreal(2)=yj
+ Rreal(3)=zj
dxj = dc_norm( 1, nres+j )
dyj = dc_norm( 2, nres+j )
dzj = dc_norm( 3, nres+j )
rij = dsqrt(rrij)
sss_ele_cut=sscale_ele(1.0d0/(rij))
sss_ele_grad=sscagrad_ele(1.0d0/(rij))
+! sss_ele_cut=1.0d0
+! sss_ele_grad=0.0d0
! print *,sss_ele_cut,sss_ele_grad,&
! 1.0d0/(rij),r_cut_ele,rlamb_ele
if (sss_ele_cut.le.0.0) cycle
sigder = fac * sigder
! fac = rij * fac
! Calculate distance derivative
- gg(1) = fac*sss_ele_cut+evdwij*sss_ele_grad
- gg(2) = fac*sss_ele_cut+evdwij*sss_ele_grad
- gg(3) = fac*sss_ele_cut+evdwij*sss_ele_grad
+ gg(1) = fac*sss_ele_cut
+ gg(2) = fac*sss_ele_cut
+ gg(3) = fac*sss_ele_cut
! if (b2.gt.0.0) then
fac = chis1 * sqom1 + chis2 * sqom2 &
- 2.0d0 * chis12 * om1 * om2 * om12
dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
dbot = 12.0d0 * b4cav * bat * Lambf
- dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow*sss_ele_cut&
- +Fcav*sss_ele_grad
- Fcav=Fcav*sss
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow*sss_ele_cut
dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
dbot = 12.0d0 * b4cav * bat * Chif
eagle = Lambf * pom
!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
gvdwx(k,i) = gvdwx(k,i) &
- - (( dFdR + gg(k) ) * pom)
+ - (( dFdR + gg(k) ) * pom)&
+ -sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
!c! & - ( dFdR * pom )
pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
gvdwx(k,j) = gvdwx(k,j) &
- + (( dFdR + gg(k) ) * pom)
+ + (( dFdR + gg(k) ) * pom) &
+ +sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
+
!c! & + ( dFdR * pom )
gvdwc(k,i) = gvdwc(k,i) &
- - (( dFdR + gg(k) ) * ertail(k))
+ - (( dFdR + gg(k) ) * ertail(k)) &
+ -sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
+
!c! & - ( dFdR * ertail(k))
gvdwc(k,j) = gvdwc(k,j) &
- + (( dFdR + gg(k) ) * ertail(k))
+ + (( dFdR + gg(k) ) * ertail(k)) &
+ +sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
+
!c! & + ( dFdR * ertail(k))
gg(k) = 0.0d0
! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
END DO
-
+
!c! Compute head-head and head-tail energies for each state
! endif
! isel=0
+! if (isel.eq.2) isel=0
+! if (isel.eq.3) isel=0
+! if (iabs(Qj).eq.1) isel=0
+! nstate(itypi,itypj)=1
IF (isel.eq.0) THEN
!c! No charges - do nothing
eheadtail = 0.0d0
CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
END IF
END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav
- evdw = evdw + Fcav + eheadtail
+ evdw = evdw + Fcav*sss_ele_cut + eheadtail*sss_ele_cut
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,&
use calc_data
use comm_momo
real (kind=8) :: facd3, facd4, federmaus, adler,&
- Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap,sgrad
! integer :: k
!c! Epol and Gpol analytical parameters
alphapol1 = alphapol(itypi,itypj)
!c! Coulomb electrostatic interaction
Ecl = (332.0d0 * Qij) / Rhead
!c! derivative of Ecl is Gcl...
- dGCLdR = (-332.0d0 * Qij ) / Rhead_sq*sss_ele_cut+ECL*sss_ele_grad
+ dGCLdR = (-332.0d0 * Qij ) / Rhead_sq*sss_ele_cut
dGCLdOM1 = 0.0d0
dGCLdOM2 = 0.0d0
dGCLdOM12 = 0.0d0
- ECL=ECL*sss_ele_grad
ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
debkap=debaykap(itypi,itypj)
-(332.0d0 * Qij *&
(dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
- dGGBdR = dGGBdFGB * dFGBdR*sss_ele_cut+Egb*sss_ele_grad
- Egb=Egb*sss_ele_cut
+ dGGBdR = dGGBdFGB * dFGBdR*sss_ele_cut
!c!-------------------------------------------------------------------
!c! Fisocav - isotropic cavity creation term
!c! or "how much energy it costs to put charged head in water"
!c! Derivative of Fisocav is GCV...
dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
dbot = 12.0d0 * al4 * pom ** 11.0d0
- dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig*sss_ele_cut
!c!-------------------------------------------------------------------
!c! Epol
!c! Polarization energy - charged heads polarize hydrophobic "neck"
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)))*sss_ele_cut+&
- (ELJ+epol)*sss_ele_grad
- epol=epol*sss_ele_cut
- Elj=Elj*sss_ele_cut
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut
!c!-------------------------------------------------------------------
!c! Return the results
!c! These things do the dRdX derivatives, that is
facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
condor = (erhead_tail(k,2) + &
facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
-
+ sgrad=(Ecl+Egb+Epol+Fisocav+Elj)*sss_ele_grad*rreal(k)*rij
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
gvdwx(k,i) = gvdwx(k,i) &
- dGCLdR * pom&
- dPOLdR1 * hawk&
- dPOLdR2 * (erhead_tail(k,2)&
-facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
- - dGLJdR * pom
+ - dGLJdR * pom-sgrad
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom&
+ dGGBdR * pom+ dGCVdR * pom&
+ dPOLdR1 * (erhead_tail(k,1)&
-facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))&
- + dPOLdR2 * condor + dGLJdR * pom
+ + dPOLdR2 * condor + dGLJdR * pom+sgrad
gvdwc(k,i) = gvdwc(k,i) &
- dGCLdR * erhead(k)&
- dGCVdR * erhead(k)&
- dPOLdR1 * erhead_tail(k,1)&
- dPOLdR2 * erhead_tail(k,2)&
- - dGLJdR * erhead(k)
+ - dGLJdR * erhead(k)-sgrad
gvdwc(k,j) = gvdwc(k,j) &
+ dGCLdR * erhead(k) &
+ dGCVdR * erhead(k) &
+ dPOLdR1 * erhead_tail(k,1) &
+ dPOLdR2 * erhead_tail(k,2)&
- + dGLJdR * erhead(k)
+ + dGLJdR * erhead(k)+sgrad
END DO
RETURN
double precision dcosom1(3),dcosom2(3)
!c! used in Epol derivatives
double precision facd3, facd4
- double precision federmaus, adler
+ double precision federmaus, adler,sgrad
integer istate,ii,jj
real (kind=8) :: Fgb
! print *,"CALLING EQUAD"
!c! Ecl = 0.0d0
!c! write (*,*) "Ecl = ", Ecl
!c! derivative of Ecl is Gcl...
- dGCLdR = (-332.0d0 * Qij ) / (Rhead_sq * eps_in)*sss_ele_cut+ECL*sss_ele_grad
- ECL=ecl*sss_ele_cut
+ dGCLdR = (-332.0d0 * Qij ) / (Rhead_sq * eps_in)
!c! dGCLdR = 0.0d0
dGCLdOM1 = 0.0d0
dGCLdOM2 = 0.0d0
dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb)
dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )&
/ ( 2.0d0 * Fgb )
- dGGBdR = dGGBdFGB * dFGBdR*sss_ele_cut+Egb*sss_ele_grad
- Egb=Egb*sss_ele_cut
+ dGGBdR = dGGBdFGB * dFGBdR
!c! dGGBdR = 0.0d0
!c!-------------------------------------------------------------------
!c! Fisocav - isotropic cavity creation term
FisoCav = top / bot
dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
dbot = 12.0d0 * al4 * pom ** 11.0d0
- dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig*sss_ele_cut+FisoCav*sss_ele_grad
- FisoCav=FisoCav*sss_ele_cut
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig
!c! dGCVdR = 0.0d0
!c!-------------------------------------------------------------------
dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
* ( 2.0d0 - 0.5d0 * ee2) ) &
/ ( 2.0d0 * fgb2 )
- dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
!c! dPOLdR1 = 0.0d0
- dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
!c! dPOLdR2 = 0.0d0
dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
!c! dPOLdOM1 = 0.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)))*sss_ele_cut+&
- (epol+Elj)*sss_ele_grad
- Elj=Elj*sss_ele_cut
- epol=epol*sss_ele_cut
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
!c! dGLJdR = 0.0d0
!c!-------------------------------------------------------------------
!c! Equad
Equad = fac * Beta1
!c! Equad = 0.0d0
!c! derivative of Equad...
- dQUADdR = ((2.5d0 * Wqd * Beta1) / (Fgb**6.0d0)) * dFGBdR*sss_ele_cut&
- + Equad*sss_ele_grad
- Equad=Equad*sss_ele_cut
+ dQUADdR = ((2.5d0 * Wqd * Beta1) / (Fgb**6.0d0)) * dFGBdR
!c! dQUADdR = 0.0d0
dQUADdOM1 = fac* (-75.0d0*om1 + 315.0d0*om1*sqom2 - 45.0d0*om2*om12)
!c! dQUADdOM1 = 0.0d0
DO k = 1, 3
dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k))
dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
- tuna(k) = eom1 * dcosom1(k) + eom2 * dcosom2(k)*sss_ele_cut
+ tuna(k) = eom1 * dcosom1(k) + eom2 * dcosom2(k)
END DO
!c! Radial stuff
DO k = 1, 3
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
!c! this acts on hydrophobic center of interaction
+! sgrad=sss_ele_grad*(Ecl+Egb+FisoCav+epol+Elj)*rij*rreal(k)
gheadtail(k,1,1) = gheadtail(k,1,1) &
- dGCLdR * pom &
- dGGBdR * pom &
DO l = 1, 4
gheadtail(k,l,2) = gheadtail(k,l,2) / eheadtail
END DO
- gvdwx(k,i) = gvdwx(k,i) + gheadtail(k,1,2)
- gvdwx(k,j) = gvdwx(k,j) + gheadtail(k,2,2)
- gvdwc(k,i) = gvdwc(k,i) + gheadtail(k,3,2)
- gvdwc(k,j) = gvdwc(k,j) + gheadtail(k,4,2)
+ gvdwx(k,i) = gvdwx(k,i) + gheadtail(k,1,2)*sss_ele_cut
+ gvdwx(k,j) = gvdwx(k,j) + gheadtail(k,2,2)*sss_ele_cut
+ gvdwc(k,i) = gvdwc(k,i) + gheadtail(k,3,2)*sss_ele_cut
+ gvdwc(k,j) = gvdwc(k,j) + gheadtail(k,4,2)*sss_ele_cut
DO l = 1, 4
gheadtail(k,l,1) = 0.0d0
gheadtail(k,l,2) = 0.0d0
END DO
END DO
eheadtail = (-dlog(eheadtail)) / betaT
+ do k=1,3
+ gvdwx(k,i) = gvdwx(k,i) - eheadtail*sss_ele_grad*rreal(k)*rij
+ gvdwx(k,j) = gvdwx(k,j) + eheadtail*sss_ele_grad*rreal(k)*rij
+ gvdwc(k,i) = gvdwc(k,i) - eheadtail*sss_ele_grad*rreal(k)*rij
+ gvdwc(k,j) = gvdwc(k,j) + eheadtail*sss_ele_grad*rreal(k)*rij
+ enddo
dPOLdOM1 = 0.0d0
dPOLdOM2 = 0.0d0
dQUADdOM1 = 0.0d0
dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
* (2.0d0 - 0.5d0 * ee1) ) &
/ (2.0d0 * fgb1)
- dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut+epol*sss_ele_grad
- epol=epol*sss_ele_cut
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut
+! epol=epol*sss_ele_cut
!c! dPOLdR1 = 0.0d0
dPOLdOM1 = 0.0d0
dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
gvdwx(k,i) = gvdwx(k,i) &
- - dPOLdR1 * hawk
+ - dPOLdR1 * hawk-epol*sss_ele_grad*rreal(k)*rij
gvdwx(k,j) = gvdwx(k,j) &
+ dPOLdR1 * (erhead_tail(k,1) &
- -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))&
+ +epol*sss_ele_grad*rreal(k)*rij
- gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1)
- gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1)
+ gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1)&
+ -epol*sss_ele_grad*rreal(k)*rij
+ gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1)&
+ +epol*sss_ele_grad*rreal(k)*rij
END DO
RETURN
dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
* (2.0d0 - 0.5d0 * ee2) ) &
/ (2.0d0 * fgb2)
- dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut+epol*sss_ele_grad
- epol=epol*sss_ele_cut
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut
+! epol=epol*sss_ele_cut
!c! dPOLdR2 = 0.0d0
dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
!c! dPOLdOM1 = 0.0d0
gvdwx(k,i) = gvdwx(k,i) &
- dPOLdR2 * (erhead_tail(k,2) &
- -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ -epol*sss_ele_grad*rreal(k)*rij
gvdwx(k,j) = gvdwx(k,j) &
- + dPOLdR2 * condor
+ + dPOLdR2 * condor+epol*sss_ele_grad*rreal(k)*rij
+
gvdwc(k,i) = gvdwc(k,i) &
- - dPOLdR2 * erhead_tail(k,2)
+ - dPOLdR2 * erhead_tail(k,2)-epol*sss_ele_grad*rreal(k)*rij
+
gvdwc(k,j) = gvdwc(k,j) &
- + dPOLdR2 * erhead_tail(k,2)
+ + dPOLdR2 * erhead_tail(k,2)+epol*sss_ele_grad*rreal(k)*rij
+
END DO
RETURN
SUBROUTINE eqd(Ecl,Elj,Epol)
use calc_data
use comm_momo
- double precision facd4, federmaus,ecl,elj,epol
+ double precision facd4, federmaus,ecl,elj,epol,sgrad
alphapol1 = alphapol(itypi,itypj)
w1 = wqdip(1,itypi,itypj)
w2 = wqdip(2,itypi,itypj)
Ecl = sparrow / Rhead**2.0d0 &
- hawk / Rhead**4.0d0
dGCLdR = (- 2.0d0 * sparrow / Rhead**3.0d0 &
- + 4.0d0 * hawk / Rhead**5.0d0)*sss_ele_cut+Ecl*sss_ele_grad
- Ecl=Ecl*sss_ele_cut
+ + 4.0d0 * hawk / Rhead**5.0d0)*sss_ele_cut
!c! dF/dom1
dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
!c! dF/dom2
dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
* (2.0d0 - 0.5d0 * ee1) ) &
/ (2.0d0 * fgb1)
- dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut+epol*sss_ele_grad
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut
!c! dPOLdR1 = 0.0d0
dPOLdOM1 = 0.0d0
dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
!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)))*sss_ele_cut+elj*sss_ele_grad
- Elj=Elj*sss_ele_cut
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut
DO k = 1, 3
erhead(k) = Rhead_distance(k)/Rhead
erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
DO k = 1, 3
hawk = (erhead_tail(k,1) + &
facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
-
+ sgrad=(epol+elj+ecl)*sss_ele_grad*rreal(k)*rij
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
gvdwx(k,i) = gvdwx(k,i) &
- dGCLdR * pom&
- dPOLdR1 * hawk &
- - dGLJdR * pom
-
+ - dGLJdR * pom &
+ -sgrad
+
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
gvdwx(k,j) = gvdwx(k,j) &
+ dGCLdR * pom &
+ dPOLdR1 * (erhead_tail(k,1) &
-facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
- + dGLJdR * pom
+ + dGLJdR * pom+sgrad
gvdwc(k,i) = gvdwc(k,i) &
- dGCLdR * erhead(k) &
- dPOLdR1 * erhead_tail(k,1) &
- - dGLJdR * erhead(k)
+ - dGLJdR * erhead(k)-sgrad
gvdwc(k,j) = gvdwc(k,j) &
+ dGCLdR * erhead(k) &
+ dPOLdR1 * erhead_tail(k,1) &
- + dGLJdR * erhead(k)
+ + dGLJdR * erhead(k)+sgrad
END DO
RETURN
use comm_momo
use calc_data
- double precision facd3, adler,ecl,elj,epol
+ double precision facd3, adler,ecl,elj,epol,sgrad
alphapol2 = alphapol(itypj,itypi)
w1 = wqdip(1,itypi,itypj)
w2 = wqdip(2,itypi,itypj)
!c! derivative of ecl is Gcl
!c! dF/dr part
dGCLdR =sss_ele_cut*(- 2.0d0 * sparrow / Rhead**3.0d0 &
- + 4.0d0 * hawk / Rhead**5.0d0)+Ecl*sss_ele_grad
+ + 4.0d0 * hawk / Rhead**5.0d0)
!c! dF/dom1
dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
!c! dF/dom2
* (2.0d0 - 0.5d0 * ee2) ) &
/ (2.0d0 * fgb2)
dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut
- epol=epol*sss_ele_cut
!c! dPOLdR2 = 0.0d0
dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
!c! dPOLdOM1 = 0.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)))*sss_ele_cut+Elj*sss_ele_grad
- elj=elj*sss_ele_cut
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut
!c!-------------------------------------------------------------------
!c! Return the results
!c! (see comments in Eqq)
DO k = 1, 3
condor = (erhead_tail(k,2) &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
-
+ sgrad=(epol+elj+ecl)*sss_ele_grad*rreal(k)*rij
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
gvdwx(k,i) = gvdwx(k,i) &
- dGCLdR * pom &
- dPOLdR2 * (erhead_tail(k,2) &
-facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
- - dGLJdR * pom
+ - dGLJdR * pom-sgrad
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
gvdwx(k,j) = gvdwx(k,j) &
+ dGCLdR * pom &
+ dPOLdR2 * condor &
- + dGLJdR * pom
+ + dGLJdR * pom+sgrad
gvdwc(k,i) = gvdwc(k,i) &
- dGCLdR * erhead(k) &
- dPOLdR2 * erhead_tail(k,2) &
- - dGLJdR * erhead(k)
+ - dGLJdR * erhead(k)-sgrad
gvdwc(k,j) = gvdwc(k,j) &
+ dGCLdR * erhead(k) &
+ dPOLdR2 * erhead_tail(k,2) &
- + dGLJdR * erhead(k)
+ + dGLJdR * erhead(k)+sgrad
END DO
RETURN
c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
* (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
- dGCLdR = (c1 - c2)*sss_ele_cut+ECL*sss_ele_grad
- ECL=ECL*sss_ele_cut
+ dGCLdR = (c1 - c2)*sss_ele_cut!+ECL*sss_ele_grad
!c! dECL/dom1
c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
DO k = 1, 3
pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
- gvdwx(k,i) = gvdwx(k,i) - dGCLdR * pom
+ gvdwx(k,i) = gvdwx(k,i)- dGCLdR * pom-(ecl*sss_ele_grad*Rreal(k)*rij)
pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
- gvdwx(k,j) = gvdwx(k,j) + dGCLdR * pom
+ gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom+(ecl*sss_ele_grad*Rreal(k)*rij)
- gvdwc(k,i) = gvdwc(k,i) - dGCLdR * erhead(k)
- gvdwc(k,j) = gvdwc(k,j) + dGCLdR * erhead(k)
+ gvdwc(k,i) = gvdwc(k,i)- dGCLdR * erhead(k)-(ecl*sss_ele_grad*Rreal(k)*rij)
+ gvdwc(k,j) = gvdwc(k,j)+ dGCLdR * erhead(k)+(ecl*sss_ele_grad*Rreal(k)*rij)
END DO
RETURN
END SUBROUTINE edd
yi=c(2,nres+i)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
- dxi=dc_norm(1,nres+i)
- dyi=dc_norm(2,nres+i)
- dzi=dc_norm(3,nres+i)
+ dxi=dc_norm(1,i)
+ dyi=dc_norm(2,i)
+ dzi=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
return
end subroutine make_cat_pep_list
+ subroutine make_lip_pep_list
+ include 'mpif.h'
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
+ real(kind=8) :: xmedj,ymedj,zmedj,sslipi,ssgradlipi,faclipij2,sslipj,ssgradlipj
+ real(kind=8) :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi
+ real(kind=8) :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj
+ real(kind=8) :: xja,yja,zja
+ integer:: contlistmartpi(300*nres),contlistmartpj(300*nres)
+ integer:: contlistmartsci(250*nres),contlistmartscj(250*nres)
+
+
+! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_martsc,&
+ ilist_martp,k,itmp
+ integer displ(0:nprocs),i_ilist_martsc(0:nprocs),ierr,&
+ i_ilist_martp(0:nprocs)
+! write(iout,*),"START make_pp",iatel_s,iatel_e,r_cut_ele+r_buff_list
+ ilist_martp=0
+ ilist_martsc=0
+
+
+ r_buff_list=6.0
+ itmp=0
+ do i=1,3
+ 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).or.(itypi.eq.10)) cycle
+! itypi1=iabs(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)
+ dxi=dc_norm(1,i)
+ dyi=dc_norm(2,i)
+ dzi=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)
+
+! dsci_inv=vbld_inv(i+nres)
+ do j=itmp+1,itmp+nres_molec(4)
+ 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)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
+ xja=boxshift(xj-xmedi,boxxsize)
+ yja=boxshift(yj-ymedi,boxysize)
+ zja=boxshift(zj-zmedi,boxzsize)
+ dist_init=xja**2+yja**2+zja**2
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_martp=ilist_martp+1
+! this can be substituted by cantor and anti-cantor
+ contlistmartpi(ilist_martp)=i
+ contlistmartpj(ilist_martp)=j
+ endif
+ xja=boxshift(xj-xi,boxxsize)
+ yja=boxshift(yj-yi,boxysize)
+ zja=boxshift(zj-zi,boxzsize)
+ dist_init=xja**2+yja**2+zja**2
+ if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
+! Here the list is created
+ ilist_martsc=ilist_martsc+1
+! this can be substituted by cantor and anti-cantor
+! write(iout,*) "have contact",i,j,ilist_martsc
+ contlistmartsci(ilist_martsc)=i
+ contlistmartscj(ilist_martsc)=j
+! write(iout,*) "have contact2",i,j,ilist_martsc,&
+! contlistmartsci(ilist_martsc),contlistmartscj(ilist_martsc)
+ endif
+! enddo
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_catsctran,ilist_catptran,&
+ ilist_catscnorm,ilist_catpnorm,ilist_catscang
+
+ do i=1,ilist_catsctran
+ write (iout,*) i,contlistcatsctrani(i),contlistcatsctranj(i),&
+ itype(j,contlistcatsctranj(i))
+ enddo
+ do i=1,ilist_catptran
+ write (iout,*) i,contlistcatptrani(i),contlistcatsctranj(i)
+ enddo
+ do i=1,ilist_catscnorm
+ write (iout,*) i,contlistcatscnormi(i),contlistcatscnormj(i)
+ enddo
+ do i=1,ilist_catpnorm
+ write (iout,*) i,contlistcatpnormi(i),contlistcatscnormj(i)
+ enddo
+ do i=1,ilist_catscang
+ write (iout,*) i,contlistcatscangi(i),contlistcatscangi(i)
+ enddo
+
+
+#endif
+ if (nfgtasks.gt.1)then
+
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ call MPI_Reduce(ilist_martsc,g_ilist_martsc,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_martsc,1,MPI_INTEGER,&
+ i_ilist_martsc,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_martsc(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistmartsci,ilist_martsc,MPI_INTEGER,&
+ newcontlistmartsci,i_ilist_martsc,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistmartscj,ilist_martsc,MPI_INTEGER,&
+ newcontlistmartscj,i_ilist_martsc,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_martsc,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistmartsci,g_ilist_martsc,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistmartscj,g_ilist_martsc,MPI_INT,king,FG_COMM,IERR)
+
+
+
+ call MPI_Reduce(ilist_martp,g_ilist_martp,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_martp,1,MPI_INTEGER,&
+ i_ilist_martp,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_martp(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistmartpi,ilist_martp,MPI_INTEGER,&
+ newcontlistmartpi,i_ilist_martp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistmartpj,ilist_martp,MPI_INTEGER,&
+ newcontlistmartpj,i_ilist_martp,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_martp,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistmartpi,g_ilist_martp,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistmartpj,g_ilist_martp,MPI_INT,king,FG_COMM,IERR)
+
+
+
+ else
+ g_ilist_martsc=ilist_martsc
+ g_ilist_martp=ilist_martp
+
+
+ do i=1,ilist_martsc
+ newcontlistmartsci(i)=contlistmartsci(i)
+ newcontlistmartscj(i)=contlistmartscj(i)
+ enddo
+ do i=1,ilist_martp
+ newcontlistmartpi(i)=contlistmartpi(i)
+ newcontlistmartpj(i)=contlistmartpj(i)
+ enddo
+ endif
+ call int_bounds(g_ilist_martsc,g_listmartsc_start,g_listmartsc_end)
+ call int_bounds(g_ilist_martp,g_listmartp_start,g_listmartp_end)
+! print *,"TUTU",g_listcatscang_start,g_listcatscang_end,i,j,g_ilist_catscangf,myrank
+
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",ilist_catsctran,ilist_catptran, &
+ ilist_catscnorm,ilist_catpnorm
+
+ do i=1,g_ilist_catsctran
+ write (iout,*) i,newcontlistcatsctrani(i),newcontlistcatsctranj(i)
+ enddo
+ do i=1,g_ilist_catptran
+ write (iout,*) i,newcontlistcatptrani(i),newcontlistcatsctranj(i)
+ enddo
+ do i=1,g_ilist_catscnorm
+ write (iout,*) i,newcontlistcatscnormi(i),newcontlistcatscnormj(i)
+ enddo
+ do i=1,g_ilist_catpnorm
+ write (iout,*) i,newcontlistcatpnormi(i),newcontlistcatscnormj(i)
+ enddo
+ do i=1,g_ilist_catscang
+ write (iout,*) i,newcontlistcatscangi(i),newcontlistcatscangj(i)
+#endif
+ return
+ end subroutine make_lip_pep_list
+
+
subroutine make_cat_cat_list
include 'mpif.h'
real(kind=8) :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
#endif
-!--------------------------------------------------------------------------
+
+!-----------LIPID-MARTINI-UNRES-PROTEIN
+
+! new for K+
+ subroutine elip_prot(evdw)
+! subroutine emart_prot2(emartion_prot)
+ use calc_data
+ use comm_momo
+
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi1,subchap,isel,itmp
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,aa,bb
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip,alpha_sco
+ integer :: ii,ki
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: escpho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,egb
+ real(kind=8) :: Fisocav,ECL,Elj,Equad,Epol,eheadtail,&
+ Lambf,&
+ Chif,ChiLambf,Fcav,dFdR,dFdOM1,&
+ emartions_prot_amber,dFdOM2,dFdL,dFdOM12,&
+ federmaus,&
+ d1i,d1j
+! real(kind=8),dimension(3,2)::erhead_tail
+! real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead,Rtail_distance
+ real(kind=8) :: facd4, adler, Fgb, facd3
+ integer troll,jj,istate
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ real(kind=8) ::locbox(3)
+ locbox(1)=boxxsize
+ locbox(2)=boxysize
+ locbox(3)=boxzsize
+
+ evdw=0.0D0
+ if (nres_molec(4).eq.0) return
+ eps_out=80.0d0
+! sss_ele_cut=1.0d0
+
+ itmp=0
+ 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
+ do ki=g_listmartsc_start,g_listmartsc_end
+ i=newcontlistmartsci(ki)
+ j=newcontlistmartscj(ki)
+
+! print *,"I am in EVDW",i
+ itypi=iabs(itype(i,1))
+
+! if (i.ne.47) 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)
+ 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=vbld_inv(i+nres)
+! do j=itmp+1,itmp+nres_molec(5)
+
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,4))
+ if ((itypj.gt.ntyp_molec(4))) cycle
+ CALL elgrad_init_mart(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+! print *,i,j,"after elgrad"
+ dscj_inv=0.0
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+
+ call to_box(xj,yj,zj)
+! write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,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)
+! write(iout,*) "xj,yj,zj", xj,yj,zj,boxxsize
+ rreal(1)=xj
+ rreal(2)=yj
+ rreal(3)=zj
+ dxj=0.0
+ dyj=0.0
+ dzj=0.0
+! dxj = dc_norm( 1, nres+j )
+! dyj = dc_norm( 2, nres+j )
+! dzj = dc_norm( 3, nres+j )
+
+ itypi = itype(i,1)
+ itypj = itype(j,4)
+! 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 = chi1mart(itypi,itypj)
+ chis1 = chis1mart(itypi,itypj)
+ chip1 = chipp1mart(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 = sigmap1mart(itypi,itypj)
+ sig2=0.0d0
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+
+! b1cav=0.0d0
+! b2cav=0.0d0
+! b3cav=0.0d0
+! b4cav=0.0d0
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabmart(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+nres)
+ ctail(k,2)=c(k,j)
+ END DO
+ call to_box(ctail(1,1),ctail(2,1),ctail(3,1))
+ call to_box(ctail(1,2),ctail(2,2),ctail(3,2))
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ do k=1,3
+ Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+ enddo
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail lomartion and distance calculations
+! dhead1
+ d1 = dheadmart(1, 1, itypi, itypj)
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! lomartion of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publimartions for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j)
+ enddo
+ call to_box(chead(1,1),chead(2,1),chead(3,1))
+ call to_box(chead(1,2),chead(2,2),chead(3,2))
+! write(iout,*) "TEST",chead(1,1),chead(2,1),chead(3,1),dc_norm(k, i+nres),d1
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ do k=1,3
+ Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k))
+ 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
+ Fisocav=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)
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
+ 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
+ if (evdw.gt.1.0d6) then
+ write (*,'(2(1x,a3,i3),7f7.2)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq
+ write(*,*) facsig,faceps1_inv,om1,chiom1,chi1
+ write(*,*) "ANISO?!",chi1
+!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+! Equad,evdwij+Fcav+eheadtail,evdw
+ endif
+
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq_mart(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_mart(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*sss_ele_cut
+!#endif
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! Calculate distance derivative
+ gg(1) = fac
+!*sss_ele_cut+evdwij*sss_ele_grad
+ gg(2) = fac
+!*sss_ele_cut+evdwij*sss_ele_grad
+ gg(3) = fac
+!*sss_ele_cut+evdwij*sss_ele_grad
+! print *,"GG(1),distance grad",gg(1)
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ 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+nres) )
+ erdxj = scalar( ertail(1), dC_norm(1,j) )
+ facd1 = dtailmart(1,itypi,itypj) * vbld_inv(i+nres)
+ facd2 = dtailmart(2,itypi,itypj) * vbld_inv(j)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ - (( dFdR + gg(k) ) * pom)*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))*sss_ele_cut&
+ +(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ gg(k) = 0.0d0
+ ENDDO
+!c! Compute head-head and head-tail energies for each state
+!! if (.false.) then ! turn off electrostatic
+ isel = iabs(Qi)+iabs(Qj)
+ if ((itype(j,4).gt.4).and.(itype(j,4).lt.14)) isel=isel+2
+! isel=0
+! if (isel.eq.2) isel=0
+ IF (isel.le.1) THEN
+ eheadtail = 0.0d0
+ ELSE IF (isel.eq.3) THEN
+ if (iabs(Qj).eq.1) then
+ CALL edq_mart(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+ else
+ 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
+ call eqd_mart(ecl,elj,epol)
+ eheadtail = ECL + elj + epol
+ endif
+ ELSE IF ((isel.eq.2)) THEN
+ if (iabs(Qi).ne.1) then
+ eheadtail=0.0d0
+ else
+ 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
+ CALL eqq_mart(Ecl,Egb,Epol,Fisocav,Elj)
+ eheadtail = ECL + Egb + Epol + Fisocav + Elj
+ endif
+ ELSE IF (isel.eq.4) then
+ call edd_mart(ecl)
+ eheadtail = ECL
+ ENDIF
+! write(iout,*) "not yet implemented",j,itype(j,5)
+!! endif ! turn off electrostatic
+ evdw = evdw + (Fcav + eheadtail)*sss_ele_cut
+! if (evdw.gt.1.0d6) then
+! write (*,'(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
+! endif
+
+ 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 (energy_dec) write(iout,*) "FCAV", &
+ sig1,sig2,b1cav,b2cav,b3cav,b4cav
+! print *,"before sc_grad_mart", i,j, gradpepmart(1,j)
+! iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad_mart
+! print *,"after sc_grad_mart", i,j, gradpepmart(1,j)
+
+! 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
+ 17 continue
+! go to 23
+! do i=ibond_start,ibond_end
+
+ do ki=g_listmartp_start,g_listmartp_end
+ i=newcontlistmartpi(ki)
+ j=newcontlistmartpj(ki)
+
+! 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
+ call to_box(xi,yi,zi)
+ 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)
+
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,4))
+ if ((itypj.gt.ntyp_molec(4))) cycle
+ CALL elgrad_init_mart_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+ dscj_inv=0.0
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ rreal(1)=xj
+ rreal(2)=yj
+ rreal(3)=zj
+
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+
+ 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,4)
+! 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 = chi1mart(itypi,itypj)
+ chis1 = chis1mart(itypi,itypj)
+ chip1 = chipp1mart(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 = sigmap1mart(itypi,itypj)
+ sig2=0.0
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabmart(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
+ call to_box(ctail(1,1),ctail(2,1),ctail(3,1))
+ call to_box(ctail(1,2),ctail(2,2),ctail(3,2))
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ do k=1,3
+ Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+ enddo
+
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail lomartion and distance calculations
+! dhead1
+ d1 = dheadmart(1, 1, itypi, itypj)
+! print *,"d1",d1
+! d1=0.0d0
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! lomartion of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publimartions 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)
+ ENDDO
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ call to_box(chead(1,1),chead(2,1),chead(3,1))
+ call to_box(chead(1,2),chead(2,2),chead(3,2))
+
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ do k=1,3
+ Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k))
+ 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 = 0.0d0 ! 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)
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
+ 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
+ om2=0.0d0
+ om12=0.0d0
+ 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
+! if (evdw.gt.1.0d6) then
+! write (*,'(2(1x,a3,i3),6f6.2)') &
+! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+! 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij
+!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+! Equad,evdwij+Fcav+eheadtail,evdw
+! endif
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq_mart(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_mart(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*sss_ele_cut
+!#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 )
+ dCAVdOM2=0.0d0
+ dCAVdOM12=0.0d0
+
+ 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 = dtailmart(1,itypi,itypj) * vbld_inv(i)
+ facd2 = dtailmart(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i))
+! gradpepmartx(k,i) = gradpepmartx(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)
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)*0.5d0
+ gradpepmart(k,i+1) = gradpepmart(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)*0.5d0
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))*sss_ele_cut&
+ +(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ gg(k) = 0.0d0
+ ENDDO
+!c! Compute head-head and head-tail energies for each state
+!c! Dipole-charge interactions
+ isel = 2+iabs(Qj)
+ if ((itype(j,4).gt.4).and.(itype(j,4).lt.14)) isel=isel+2
+! if (isel.eq.4) isel=0
+ if (isel.le.2) then
+ eheadtail=0.0d0
+ ELSE if (isel.eq.3) then
+ CALL edq_mart_pep(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! print *,"i,",i,eheadtail
+! eheadtail = 0.0d0
+ else
+!HERE WATER and other types of molecules solvents will be added
+! write(iout,*) "not yet implemented"
+ CALL edd_mart_pep(ecl)
+ eheadtail=ecl
+! CALL edd_mart_pep
+! eheadtail=0.0d0
+ endif
+ evdw = evdw +( Fcav + eheadtail)*sss_ele_cut
+! if (evdw.gt.1.0d6) then
+! write (*,'(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
+! endif
+ 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_mart_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
+ 23 continue
+! print *,"before leave sc_grad_mart", i,j, gradpepmart(1,nres-1)
+
+ return
+ end subroutine elip_prot
+
+ SUBROUTINE eqq_mart(Ecl,Egb,Epol,Fisocav,Elj)
+ use calc_data
+ use comm_momo
+ real (kind=8) :: facd3, facd4, federmaus, adler,&
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
+! integer :: k
+!c! Epol and Gpol analytical parameters
+ alphapol1 = alphapolmart(itypi,itypj)
+ alphapol2 = alphapolmart2(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphisomart(1,itypi,itypj)
+ al2 = alphisomart(2,itypi,itypj)
+ al3 = alphisomart(3,itypi,itypj)
+ al4 = alphisomart(4,itypi,itypj)
+ csig = (1.0d0 &
+ / dsqrt(sigiso1mart(itypi, itypj)**2.0d0 &
+ + sigiso2mart(itypi,itypj)**2.0d0))
+!c!
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(itypi,itypj)
+ Rhead_sq = Rhead * Rhead
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R1 = 0.0d0
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances needed by Epol
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+ 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! Coulomb electrostatic interaction
+ Ecl = (332.0d0 * Qij) / Rhead
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij ) / Rhead_sq
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ debkap=debaykapmart(itypi,itypj)
+ if (energy_dec) write(iout,*) "egb",Qij,debkap,Fgb,a12sq,ee0
+ Egb = -(332.0d0 * Qij *&
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
+! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
+!c! Derivative of Egb is Ggb...
+ dGGBdFGB = -(-332.0d0 * Qij * &
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+ -(332.0d0 * Qij *&
+ (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
+ dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
+ dGGBdR = dGGBdFGB * dFGBdR
+!c!-------------------------------------------------------------------
+!c! Fisocav - isotropic cavity creation term
+!c! or "how much energy it costs to put charged head in water"
+ pom = Rhead * csig
+ top = al1 * (dsqrt(pom) + al2 * pom - al3)
+ bot = (1.0d0 + al4 * pom**12.0d0)
+ botsq = bot * bot
+ FisoCav = top / bot
+! write (*,*) "Rhead = ",Rhead
+! write (*,*) "csig = ",csig
+! write (*,*) "pom = ",pom
+! write (*,*) "al1 = ",al1
+! write (*,*) "al2 = ",al2
+! write (*,*) "al3 = ",al3
+! write (*,*) "al4 = ",al4
+! write (*,*) "top = ",top
+! write (*,*) "bot = ",bot
+!c! Derivative of Fisocav is GCV...
+ dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
+ dbot = 12.0d0 * al4 * pom ** 11.0d0
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig
+!c!-------------------------------------------------------------------
+!c! Epol
+!c! Polarization energy - charged heads polarize hydrophobic "neck"
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR1 = ( R1 * R1 ) / MomoFac1
+ RR2 = ( R2 * R2 ) / MomoFac2
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ ee2 = exp(-( RR2 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1 )
+ fgb2 = sqrt( RR2 + a12sq * ee2 )
+ epol = 332.0d0 * eps_inout_fac * ( &
+ (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 ))
+!c! epol = 0.0d0
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)&
+ / (fgb1 ** 5.0d0)
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)&
+ / (fgb2 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1)* ( 2.0d0 - (0.5d0 * ee1) ) )&
+ / ( 2.0d0 * fgb1 )
+ dFGBdR2 = ( (R2 / MomoFac2)* ( 2.0d0 - (0.5d0 * ee2) ) )&
+ / ( 2.0d0 * fgb2 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1))&
+ * ( 2.0d0 - 0.5d0 * ee1) ) / ( 2.0d0 * fgb1 )
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2))&
+ * ( 2.0d0 - 0.5d0 * ee2) ) / ( 2.0d0 * fgb2 )
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1!*sss_ele_cut+epol*sss_ele_grad
+!c! dPOLdR1 = 0.0d0
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2!*sss_ele_cut+epol*sss_ele_grad
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+! epol=epol*sss_ele_cut
+!c! dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+!c! Lennard-Jones 6-12 interaction between heads
+ 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! These things do the dRdX derivatives, that is
+!c! allow us to change what we see from function that changes with
+!c! distance to function that changes with LOCATION (of the interaction
+!c! site)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ 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) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,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+nres)
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailmart(1,itypi,itypj) * vbld_inv(i+nres)
+ facd4 = dtailmart(2,itypi,itypj) * vbld_inv(j)
+
+!c! Now we add appropriate partial derivatives (one in each dimension)
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+ 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))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ +sss_ele_cut*(- dGCLdR * pom&
+ - dGGBdR * pom&
+ - dGCVdR * pom&
+ - dPOLdR1 * hawk&
+ - dPOLdR2 * (erhead_tail(k,2)&
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ - dGLJdR * pom)-&
+ sss_ele_grad*rij*rreal(k)*(Ecl+Egb+Epol+Fisocav+Elj)
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepmartx(k,j) = gradpepmartx(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
+
+ gradpepmart(k,i) = gradpepmart(k,i) + &
+ sss_ele_cut*(- dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k))&
+ - sss_ele_grad*rij*rreal(k)*(Ecl+Egb+Epol+Fisocav+Elj)
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) + &
+ sss_ele_cut*( dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2)&
+ + dGLJdR * erhead(k))&
+ +sss_ele_grad*rij*rreal(k)*(Ecl+Egb+Epol+Fisocav+Elj)
+ END DO
+ RETURN
+ END SUBROUTINE eqq_mart
+
+ SUBROUTINE eqd_mart(Ecl,Elj,Epol)
+ use calc_data
+ use comm_momo
+ double precision facd4, federmaus,ecl,elj,epol
+ alphapol1 = alphapolmart(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(itypi,itypj)
+! eps_head=0.0d0
+! w2=0.0d0
+! alphapol1=0.0d0
+!c!-------------------------------------------------------------------
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!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 * Qi * om1
+ hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ Ecl = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+ dGCLdR =sss_ele_cut*(-2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0)
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = 0.0d0 !
+
+!(2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+!c! epol = 0.0d0
+!c!------------------------------------------------------------------
+!c! derivative of Epol is Gpol...
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) &
+ / (fgb1 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM2 = 0.0d0 ! as om2 is 0
+! (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+! * (2.0d0 - 0.5d0 * ee1) ) &
+! / (2.0d0 * fgb1)
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut
+!c! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+! dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ 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*sss_ele_cut &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ - dGCLdR * pom&
+ - dPOLdR1 * hawk &
+ - dGLJdR * pom&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepmartx(k,j) = gradpepmartx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR1 * (erhead_tail(k,1) &
+! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+! + dGLJdR * pom
+
+
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1) &
+ - dGLJdR * erhead(k)&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dGLJdR * erhead(k)&
+ +(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ END DO
+ RETURN
+ END SUBROUTINE eqd_mart
+
+ SUBROUTINE edq_mart(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolmart(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(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
+! write(iout,*) "KURWA2",Rhead
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+ 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)*sss_ele_cut
+!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*sss_ele_cut
+!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)))*sss_ele_cut
+!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+nres)
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailmart(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))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+ - dGLJdR * pom&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepmartx(k,j) = gradpepmartx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k)&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)&
+ +(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_mart
+
+ SUBROUTINE edq_mart_pep(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolmart(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(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 *,"CO2", itypi,itypj
+! 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)*sss_ele_cut
+!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*sss_ele_cut
+!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*sss_ele_cut &
+ * (((-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) )
+ facd1 = d1 * vbld_inv(i+1)
+ DO k = 1, 3
+ pom = facd1*(erhead(k)-erdxi*dC_norm(k,i))
+! gradpepmartx(k,i) = gradpepmartx(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))
+! gradpepmartx(k,j) = gradpepmartx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+ gradpepmart(k,i) = gradpepmart(k,i)+pom*(dGCLdR+dGLJdR)
+ gradpepmart(k,i+1) = gradpepmart(k,i+1)-pom*(dGCLdR+dGLJdR)
+
+ gradpepmart(k,i) = gradpepmart(k,i) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+ gradpepmart(k,i+1) = gradpepmart(k,i+1) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)&
+ +(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_mart_pep
+!--------------------------------------------------------------------------
+
+ SUBROUTINE edd_mart(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+! w2=0.0d0
+!c!-------------------------------------------------------------------
+!c! ECL
+! print *,"om1",om1,om2,om12
+ fac = - 3.0d0 * om1 !after integer and simplify
+ c1 = (w1 / (Rhead**3.0d0)) * fac
+ c2 = (w2 / Rhead ** 6.0d0) &
+ * (4.0d0 + 6.0d0*sqom1 ) !after integration and simplifimartion
+ ECL = c1 - c2
+!c! dervative of ECL is GCL...
+!c! dECL/dr
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + 6.0d0*sqom1)
+ dGCLdR = (c1 - c2)*sss_ele_cut
+!c! dECL/dom1
+ c1 = (-3.0d0 * w1) / (Rhead**3.0d0)
+ c2 = (12.0d0 * w2*om1) / (Rhead**6.0d0)
+ dGCLdOM1 = c1 - c2
+!c! dECL/dom2
+! c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
+ c1=0.0 ! this is because om2 is 0
+! c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+! * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
+ c2=0.0 !om is 0
+ dGCLdOM2 = c1 - c2
+!c! dECL/dom12
+! c1 = w1 / (Rhead ** 3.0d0)
+ c1=0.0d0 ! this is because om12 is 0
+! c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ c2=0.0d0 !om12 is 0
+ dGCLdOM12 = c1 - c2
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k= 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ DO k = 1, 3
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gradpepmartx(k,i) = gradpepmartx(k,i) - dGCLdR * pom&
+ -ecl*sss_ele_grad*rij*rreal(k)
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepmartx(k,j) = gradpepmartx(k,j) + dGCLdR * pom
+
+ gradpepmart(k,i) = gradpepmart(k,i) - dGCLdR * erhead(k)&
+ -ecl*sss_ele_grad*rij*rreal(k)
+
+ gradpepmart(k,j) = gradpepmart(k,j) + dGCLdR * erhead(k)&
+ +ecl*sss_ele_grad*rij*rreal(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edd_mart
+ SUBROUTINE edd_mart_pep(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! ECL
+ fac = (om12 - 3.0d0 * om1 * om2)
+ c1 = (w1 / (Rhead**3.0d0)) * fac
+ c2 = (w2 / Rhead ** 6.0d0) &
+ * (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+ ECL = c1 - c2
+!c! dECL/dr
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
+ dGCLdR = (c1 - c2)*sss_ele_cut
+!c! dECL/dom1
+ c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 )
+ dGCLdOM1 = c1 - c2
+!c! dECL/dom2
+ c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
+ dGCLdOM2 = c1 - c2
+ dGCLdOM2=0.0d0 ! this is because om2=0
+!c! dECL/dom12
+ c1 = w1 / (Rhead ** 3.0d0)
+ c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ dGCLdOM12 = c1 - c2
+ dGCLdOM12=0.0d0 !this is because om12=0.0
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k= 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ facd1 = d1 * vbld_inv(i)
+ facd2 = d2 * vbld_inv(j+nres)
+ DO k = 1, 3
+
+ pom = facd1*(erhead(k)-erdxi*dC_norm(k,i))
+ gradpepmart(k,i) = gradpepmart(k,i) + dGCLdR * pom
+ gradpepmart(k,i+1) = gradpepmart(k,i+1) - dGCLdR * pom
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepmartx(k,j) = gradpepmartx(k,j) + dGCLdR * pom
+
+ gradpepmart(k,i) = gradpepmart(k,i) - dGCLdR * erhead(k)*0.5d0&
+ -ECL*sss_ele_grad*rreal(k)*rij
+ gradpepmart(k,i+1) = gradpepmart(k,i+1)- dGCLdR * erhead(k)*0.5d0&
+ -ECL*sss_ele_grad*rreal(k)*rij
+
+ gradpepmart(k,j) = gradpepmart(k,j) + dGCLdR * erhead(k)&
+ +ECL*sss_ele_grad*rreal(k)*rij
+
+ END DO
+ RETURN
+ END SUBROUTINE edd_mart_pep
+
+ SUBROUTINE elgrad_init_mart(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 = itype(i,1)
+ itypj = itype(j,4)
+! print *,"in elegrad",i,j,itypi,itypj
+!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 = sigmamart( itypi,itypj )
+ chi1 = chi1mart( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1mart( 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 )
+! print *,"before dheadmart"
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1mart(itypi,itypj) * rborn2mart(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+! print *,"after dheadmart"
+ Qi = icharge(itypi)
+ Qj = ichargelipid(itypj)
+ Qij = Qi * Qj
+! print *,"after icharge"
+
+!c! chis1,2,12
+ chis1 = chis1mart(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1mart(itypi,itypj)
+ sig2 = sigmap2mart(itypi,itypj)
+! print *,"before alphasurmart"
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+ wqd = wquadmart(itypi, itypj)
+! print *,"after alphasurmar n wquad"
+!c! used by Fgb
+ eps_in = epsintabmart(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail lomartion and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtailmart(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j)!-dtailmart(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 lomartion and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! lomartion of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publimartions for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+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)
+ 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_mart
+
+ SUBROUTINE elgrad_init_mart_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,4)
+!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 = sigmamart( itypi,itypj )
+ chi1 = chi1mart( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1mart( 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 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1mart(itypi,itypj) * rborn2mart(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = 0
+ Qj = ichargelipid(itypj)
+! Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis1mart(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1mart(itypi,itypj)
+ sig2 = sigmap2mart(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+ wqd = wquadmart(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabmart(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail lomartion and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0-dtailmart(1,itypi,itypj)*dc_norm(k,i)
+ ctail(k,2)=c(k,j)!-dtailmart(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 lomartion and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! lomartion of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publimartions 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_mart_pep
+
+ subroutine sc_grad_mart
+ 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,nres+i)-om1*erij(k))
+ dcosom2(k)=rij*(dc_norm(k,j)-om2*erij(k))
+ enddo
+ do k=1,3
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))
+! print *,'gg',k,gg(k)
+ enddo
+! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut
+! write (iout,*) "gg",(gg(k),k=1,3)
+ do k=1,3
+ gradpepmartx(k,i)=gradpepmartx(k,i)-gg(k)*sss_ele_cut &
+ +(eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss_ele_cut
+
+! gradpepcatx(k,j)=gradpepcatx(k,j)+gg(k) &
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv
+
+! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!
+! Calculate the components of the gradient in DC and X
+!
+ do l=1,3
+ gradpepmart(l,i)=gradpepmart(l,i)-gg(l)*sss_ele_cut
+ gradpepmart(l,j)=gradpepmart(l,j)+gg(l)*sss_ele_cut
+ enddo
+ end subroutine sc_grad_mart
+
+ subroutine sc_grad_mart_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
+! write (iout,*) "gg",(gg(k),k=1,3)
+
+ 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)
+ gradpepmart(k,i)= gradpepmart(k,i) +sss_ele_cut*(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)
+ gradpepmart(k,i+1)= gradpepmart(k,i+1) +sss_ele_cut*(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)
+ gradpepmart(k,j)=gradpepmart(k,j)+gg(k)*sss_ele_cut
+ enddo
+ end subroutine sc_grad_mart_pep
end module energy
projects / unres4.git / blobdiff