X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=bb7d08d5afe704eb0f2bd1c4bc4afffc7d5f22aa;hb=refs%2Fheads%2FUCGM;hp=d81310f1133ab0161691d4cc17a3baa288a5ab79;hpb=91b9f78d94b96277537615722323ebe03cc0a014;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index d81310f..bb7d08d 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -142,7 +142,7 @@ ! 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) @@ -259,7 +259,7 @@ ! 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 @@ -274,8 +274,10 @@ integer ishield_listbuf(-1:nres_molec(1)), & shield_listbuf(maxcontsshi,-1:nres_molec(1)),k,j,i,iii,impishi,mojint,jjj + integer :: imatupdate2 ! print *,"I START ENERGY" imatupdate=100 + imatupdate2=100 ! if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list ! real(kind=8), dimension(:),allocatable:: fac_shieldbuf ! real(kind=8), dimension(:,:,:),allocatable:: & @@ -345,7 +347,8 @@ weights_(49)=wpeppho weights_(50)=wcatnucl weights_(56)=wcat_tran - + weights_(58)=wlip_prot + weights_(52)=wmartini ! wcatcat= weights(41) ! wcatprot=weights(42) @@ -395,8 +398,9 @@ 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) @@ -409,7 +413,7 @@ time_Bcastw=time_Bcastw+MPI_Wtime()-time00 ! call chainbuild_cart endif - print *,"itime_mat",itime_mat,imatupdate +! print *,"itime_mat",itime_mat,imatupdate if (nfgtasks.gt.1) then call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR) endif @@ -418,15 +422,26 @@ ! 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 ! print *,'Processor',myrank,' calling etotal ipot=',ipot call make_cat_pep_list +! call make_cat_cat_list endif endif endif + if (nres_molec(5).gt.0) then + if (mod(itime_mat,imatupdate2).eq.0) then +! print *, "before cat cat" +! print *,'Processor',myrank,' calling etotal ipot=',ipot +! call make_cat_pep_list + call make_cat_cat_list + endif + endif ! write (iout,*) "after make_pp_inter_list" ! print *,'Processor',myrank,' calling etotal ipot=',ipot @@ -893,7 +908,7 @@ Eafmforce=0.0d0 endif endif - print *,"before tubemode",tubemode +! print *,"before tubemode",tubemode if (tubemode.eq.1) then call calctube(etube) else if (tubemode.eq.2) then @@ -903,7 +918,7 @@ else etube=0.0d0 endif - print *, "TU JEST PRZED EHPB" +! print *, "TU JEST PRZED EHPB" call edis(ehpb) !-------------------------------------------------------- @@ -949,13 +964,13 @@ else ecation_protang=0.0d0 endif - if (nfgtasks.gt.1) then - if (fg_rank.eq.0) then +! if (nfgtasks.gt.1) then +! if (fg_rank.eq.0) then if (nres_molec(5).gt.1) call ecatcat(ecationcation) - endif - else - if (nres_molec(5).gt.1) call ecatcat(ecationcation) - endif +! endif +! else +! if (nres_molec(5).gt.1) call ecatcat(ecationcation) +! endif if (oldion.gt.0) then if (g_ilist_catpnorm.gt.0) call ecat_prot(ecation_prot) else @@ -992,6 +1007,11 @@ 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 @@ -1079,6 +1099,7 @@ 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" @@ -1128,7 +1149,7 @@ 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 @@ -1218,6 +1239,7 @@ elipelec=energia(55) ecat_prottran=energia(56) ecation_protang=energia(57) + elipidprot=energia(58) ! ecations_prot_amber=energia(50) ! energia(41)=ecation_prot @@ -1238,7 +1260,9 @@ +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 @@ -1258,7 +1282,9 @@ +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 @@ -1399,7 +1425,7 @@ 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) @@ -1459,7 +1485,7 @@ 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,& @@ -1479,7 +1505,7 @@ 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)'/ & @@ -1534,6 +1560,7 @@ '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,& @@ -1550,7 +1577,7 @@ 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)'/ & @@ -1602,6 +1629,7 @@ '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 @@ -12096,8 +12124,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) + @@ -12136,8 +12166,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) @@ -12409,7 +12440,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! +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 @@ -12720,9 +12753,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut ! write (iout,*) "gg",(gg(k),k=1,3) do k=1,3 - gradpepcatx(k,i)=gradpepcatx(k,i)-gg(k) & + gradpepcatx(k,i)=gradpepcatx(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 + +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)) & @@ -12737,8 +12770,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! Calculate the components of the gradient in DC and X ! do l=1,3 - gradpepcat(l,i)=gradpepcat(l,i)-gg(l) - gradpepcat(l,j)=gradpepcat(l,j)+gg(l) + gradpepcat(l,i)=gradpepcat(l,i)-gg(l)*sss_ele_cut + gradpepcat(l,j)=gradpepcat(l,j)+gg(l)*sss_ele_cut enddo end subroutine sc_grad_cat @@ -12758,20 +12791,21 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! 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) - gvdwc_pepbase(k,i)= gvdwc_pepbase(k,i) +0.5*(- gg(k)) & + gradpepcat(k,i)= gradpepcat(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 - gvdwc_pepbase(k,i+1)= gvdwc_pepbase(k,i+1) +0.5*(- gg(k)) & + - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0) + gradpepcat(k,i+1)= gradpepcat(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 - gradpepcat(k,j)=gradpepcat(k,j)+gg(k) + + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0) + gradpepcat(k,j)=gradpepcat(k,j)+gg(k)*sss_ele_cut enddo end subroutine sc_grad_cat_pep @@ -13420,8 +13454,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! include 'COMMON.VAR' ! include 'COMMON.CONTACTS' use comm_srutu +!#ifdef LBFGS +! use minimm, only: funcgrad +!#endif !el integer :: icall !el common /srutu/ icall +! real(kind=8) :: funcgrad real(kind=8),dimension(6) :: ggg real(kind=8),dimension(3) :: cc,xx,ddc,ddx real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres) @@ -13431,7 +13469,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! real(kind=8) :: urparm(1) !EL external fdum integer :: nf,i,j,k - real(kind=8) :: aincr,etot,etot1 + real(kind=8) :: aincr,etot,etot1,ff icg=1 nf=0 nfl=0 @@ -13443,8 +13481,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! call geom_to_var(nvar,x) call etotal(energia) etot=energia(0) +#ifdef LBFGS + ff=funcgrad(x,g) +#else !el call enerprint(energia) call gradient(nvar,x,nf,g,uiparm,urparm,fdum) +#endif icall =1 do i=1,nres write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) @@ -13537,7 +13579,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=1.0D-5 + aincr=graddelta write(iout,*) 'Calling CHECK_ECARTINT.,kupa' nf=0 icall=0 @@ -13997,8 +14039,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! include 'COMMON.VAR' ! include 'COMMON.GEO' use comm_srutu +!#ifdef LBFGS +! use minimm, only : funcgrad +!#endif !el integer :: icall !el common /srutu/ icall +! real(kind=8) :: funcgrad real(kind=8),dimension(6*nres) :: x,gana,gg !(maxvar) (maxvar=6*maxres) integer :: uiparm(1) real(kind=8) :: urparm(1) @@ -14006,7 +14052,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! character(len=6) :: key !EL external fdum integer :: i,ii,nf - real(kind=8) :: xi,aincr,etot,etot1,etot2 + real(kind=8) :: xi,aincr,etot,etot1,etot2,ff call zerograd aincr=1.0D-7 print '(a)','Calling CHECK_INT.' @@ -14032,9 +14078,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! #endif nf=1 nfl=3 +#ifdef LBFGS + ff=funcgrad(x,gana) +#else + !d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar) call gradient(nvar,x,nf,gana,uiparm,urparm,fdum) !d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp +#endif icall=1 do i=1,nvar xi=x(i) @@ -15763,7 +15814,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! #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 @@ -17713,6 +17764,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !----------------------------------------------------------------------------- ! gradient_p.F !----------------------------------------------------------------------------- +#ifndef LBFGS subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm) use io_base, only:intout,briefout @@ -17818,6 +17870,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !d write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n) return end subroutine gradient +#endif !----------------------------------------------------------------------------- subroutine func(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F @@ -18190,6 +18243,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -21763,6 +21818,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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)) @@ -21964,6 +22021,11 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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)) @@ -21972,7 +22034,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! allocate(newcontlistcatscnormj(300*nres)) allocate(newcontlistcatpnormi(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)) @@ -23608,7 +23672,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! #endif subroutine ecatcat(ecationcation) use MD_data, only: t_bath - integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj,irdiff + integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj,irdiff,& + ii real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,& r7,r4,ecationcation,k0,rcal,aa,bb,sslipi,ssgradlipi,sslipj,ssgradlipj real(kind=8) :: xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, & @@ -23626,12 +23691,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! k0 = 332.0*(2.0*2.0)/80.0 itmp=0 - do i=1,4 - itmp=itmp+nres_molec(i) - enddo -! write(iout,*) "itmp",itmp - do i=itmp+1,itmp+nres_molec(5)-1 - +! do i=1,4 +! itmp=itmp+nres_molec(i) +! enddo +! write(iout,*) "itmp",g_listcatcatnorm_start, g_listcatcatnorm_end +! do i=itmp+1,itmp+nres_molec(5)-1 + do ii=g_listcatcatnorm_start, g_listcatcatnorm_end + i=newcontlistcatcatnormi(ii) + j=newcontlistcatcatnormj(ii) + xi=c(1,i) yi=c(2,i) zi=c(3,i) @@ -23639,7 +23707,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! itypi=itype(i,5) call to_box(xi,yi,zi) call lipid_layer(xi,yi,zi,sslipi,ssgradlipi) - do j=i+1,itmp+nres_molec(5) +! do j=i+1,itmp+nres_molec(5) itypj=itype(j,5) ! print *,i,j,itypi,itypj k0 = 332.0*(ichargecat(itypi)*ichargecat(itypj))/80.0 @@ -23659,7 +23727,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! rcal =xj**2+yj**2+zj**2 ract=sqrt(rcal) if ((itypi.gt.1).or.(itypj.gt.1)) then - + if (sss2min2.eq.0.0d0) cycle + sss2min2=sscale2(ract,12.0d0,1.0d0) + sss2mingrad2=sscagrad2(ract,12.0d0,1.0d0) ! rcat0=3.472 ! epscalc=0.05 ! r06 = rcat0**6 @@ -23680,13 +23750,13 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! enddo do k=1,3 gg(k) = dEvan1Cmcat(k)+dEvan2Cmcat(k)+dEeleccat(k) - gradcatcat(k,i)=gradcatcat(k,i)-gg(k) - gradcatcat(k,j)=gradcatcat(k,j)+gg(k) + gradcatcat(k,i)=gradcatcat(k,i)-(gg(k)*sss2min2+(Evan1cat+Evan2cat+Eeleccat)*sss2mingrad2) + gradcatcat(k,j)=gradcatcat(k,j)+gg(k)*sss2min2+(Evan1cat+Evan2cat+Eeleccat)*sss2mingrad2 enddo if (energy_dec) write (iout,*) "ecatcat",i,j,Evan1cat,Evan2cat,Eeleccat,& r012,rcal**6,ichargecat(itypi)*ichargecat(itypj) ! write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj - ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat + ecationcation=ecationcation+(Evan1cat+Evan2cat+Eeleccat)*sss2min2 else !this is water part and other non standard molecules sss2min2=sscale2(ract,10.0d0,1.0d0)! cutoff for water interaction is 15A @@ -23715,10 +23785,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! gradcatcat(k,i)=gradcatcat(k,i)-gg(k)*sss2min2-sss2mingrad2*ewater*r(k)/ract gradcatcat(k,j)=gradcatcat(k,j)+gg(k)*sss2min2+sss2mingrad2*ewater*r(k)/ract enddo - if (energy_dec) write(iout,'(2f10.7,f15.7,2i5)') rdiff,ract,ecationcation,i,j + if (energy_dec) write(iout,'(2f8.2,f10.2,2i5)') rdiff,ract,ecationcation,i,j endif ! end water enddo - enddo +! enddo return end subroutine ecatcat !--------------------------------------------------------------------------- @@ -23736,7 +23806,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 + integer :: ii,ki real(kind=8) :: fracinbuf real (kind=8) :: escpho real (kind=8),dimension(4):: ener @@ -23768,7 +23838,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 i=ibond_start,ibond_end + do ki=g_listcatscnorm_start,g_listcatscnorm_end + i=newcontlistcatscnormi(ki) + j=newcontlistcatscnormj(ki) ! print *,"I am in EVDW",i itypi=iabs(itype(i,1)) @@ -23785,7 +23858,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) +! do j=itmp+1,itmp+nres_molec(5) ! Calculate SC interaction energy. itypj=iabs(itype(j,5)) @@ -23810,6 +23883,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! zj=boxshift(zj-zi,boxzsize) ! write(iout,*) "xj,yj,zj", xj,yj,zj,boxxsize + 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 ) @@ -23919,6 +23995,11 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! 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 @@ -23969,15 +24050,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! END IF !#else evdw = evdw & - + evdwij + + 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 + 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 @@ -23996,8 +24077,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf)) dbot = 12.0d0 * b4cav * bat * Lambf - dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow - + dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow*sss_ele_cut+& + Fcav*sss_ele_grad + Fcav=Fcav*sss_ele_cut dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif)) dbot = 12.0d0 * b4cav * bat * Chif eagle = Lambf * pom @@ -24037,61 +24119,47 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj) ! print *,i,itype(i,1),isel IF (isel.eq.0) THEN -!c! No charges - do nothing eheadtail = 0.0d0 - ELSE IF (isel.eq.1) THEN -!c! Nonpolar-charge interactions if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then Qi=Qi*2 Qij=Qij*2 endif - CALL enq_cat(epol) eheadtail = epol -! eheadtail = 0.0d0 - ELSE IF (isel.eq.3) THEN -!c! Dipole-charge interactions if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then Qi=Qi*2 Qij=Qij*2 endif -! write(iout,*) "KURWA0",d1 - CALL edq_cat(ecl, elj, epol) eheadtail = ECL + elj + epol -! eheadtail = 0.0d0 - ELSE IF ((isel.eq.2)) THEN - -!c! Same charge-charge interaction ( +/+ or -/- ) 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_cat(Ecl,Egb,Epol,Fisocav,Elj) eheadtail = ECL + Egb + Epol + Fisocav + Elj -! eheadtail = 0.0d0 - -! ELSE IF ((isel.eq.2.and. & -! iabs(Qi).eq.1).and. & -! nstate(itypi,itypj).ne.1) THEN -!c! Different charge-charge interaction ( +/- or -/+ ) -! if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then -! Qi=Qi*2 -! Qij=Qij*2 -! endif -! if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then -! Qj=Qj*2 -! Qij=Qij*2 -! endif -! -! CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav - else - write(iout,*) "not yet implemented",j,itype(j,5) + else ! here is water and other molecules + isel = iabs(Qi)+2 +! isel=2 +! if (isel.eq.4) isel=2 + if (isel.eq.2) then + eheadtail = 0.0d0 + else if (isel.eq.3) then + 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_cat(ecl,elj,epol) + eheadtail = ECL + elj + epol + else if (isel.eq.4) then + call edd_cat(ecl) + eheadtail = ECL + endif +! write(iout,*) "not yet implemented",j,itype(j,5) endif !! endif ! turn off electrostatic evdw = evdw + Fcav + eheadtail @@ -24118,14 +24186,18 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !c!------------------------------------------------------------------- !c! NAPISY KONCOWE END DO ! j - END DO ! i +! 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 i=ibond_start,ibond_end + + do ki=g_listcatpnorm_start,g_listcatpnorm_end + i=newcontlistcatpnormi(ki) + j=newcontlistcatpnormj(ki) ! print *,"I am in EVDW",i itypi=10 ! the peptide group parameters are for glicine @@ -24141,7 +24213,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) +! do j=itmp+1,itmp+nres_molec(5) ! Calculate SC interaction energy. itypj=iabs(itype(j,5)) @@ -24266,15 +24338,22 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dCAVdOM1 = 0.0d0 dCAVdOM2 = 0.0d0 dCAVdOM12 = 0.0d0 - dscj_inv = vbld_inv(j+nres) + 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 @@ -24318,15 +24397,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! END IF !#else evdw = evdw & - + evdwij + + 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 + 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 fac = chis1 * sqom1 + chis2 * sqom2 & - 2.0d0 * chis12 * om1 * om2 * om12 @@ -24347,20 +24426,24 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf)) dbot = 12.0d0 * b4cav * bat * Lambf - dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow - + dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow*sss_ele_cut+& + Fcav*sss_ele_grad + Fcav=Fcav*sss_ele_cut 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 = dFdL * ( dFdOM2 ) +! dCAVdOM12 = dFdL * ( dFdOM12 ) + dCAVdOM2=0.0d0 + dCAVdOM12=0.0d0 DO k= 1, 3 ertail(k) = Rtail_distance(k)/Rtail @@ -24395,8 +24478,11 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! eheadtail = 0.0d0 else !HERE WATER and other types of molecules solvents will be added - write(iout,*) "not yet implemented" +! write(iout,*) "not yet implemented" + CALL edd_cat_pep(ecl) + eheadtail=ecl ! CALL edd_cat_pep +! eheadtail=0.0d0 endif evdw = evdw + Fcav + eheadtail ! if (evdw.gt.1.0d6) then @@ -24417,7 +24503,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !c!------------------------------------------------------------------- !c! NAPISY KONCOWE END DO ! j - END DO ! i +! END DO ! i !c write (iout,*) "Number of loop steps in EGB:",ind !c energy_dec=.false. ! print *,"EVDW KURW",evdw,nres @@ -24467,6 +24553,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization do i=ibond_start,ibond_end ! cycle + if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms xi=0.5d0*(c(1,i)+c(1,i+1)) yi=0.5d0*(c(2,i)+c(2,i+1)) @@ -26780,7 +26867,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 + integer :: ii,icont real(kind=8) :: fracinbuf real (kind=8) :: escpho real (kind=8),dimension(4):: ener @@ -26801,8 +26888,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! sss_ele_cut=1.0d0 countss=0 ! print *,"EVDW KURW",evdw,nres - do i=iatsc_s,iatsc_e +! do i=iatsc_s,iatsc_e ! print *,"I am in EVDW",i + do icont=g_listscsc_start,g_listscsc_end + i=newcontlisti(icont) + j=newcontlistj(icont) + itypi=iabs(itype(i,1)) ! if (i.ne.47) cycle if (itypi.eq.ntyp1) cycle @@ -26824,8 +26915,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! ! Calculate SC interaction energy. ! - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +! do iint=1,nint_gr(i) +! do j=istart(i,iint),iend(i,iint) ! print *,"JA PIER",i,j,iint,istart(i,iint),iend(i,iint) IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN call dyn_ssbond_ene(i,j,evdwij,countss) @@ -26877,6 +26968,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 ) @@ -27003,6 +27097,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! 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)) +! 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 + !---------------------------- CALL sc_angular ! this should be in elgrad_init but om's are calculated by sc_angular @@ -27044,7 +27146,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! END IF !#else evdw = evdw & - + evdwij + + evdwij*sss_ele_cut !#endif c1 = c1 * eps1 * eps2rt**2 * eps3rt**2 @@ -27052,9 +27154,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! sigder = fac * sigder ! fac = rij * fac ! Calculate distance derivative - gg(1) = fac - gg(2) = fac - gg(3) = fac + 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 @@ -27079,8 +27181,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf)) dbot = 12.0d0 * b4cav * bat * Lambf - dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow - + 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 @@ -27107,26 +27208,33 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !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 @@ -27142,6 +27250,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! 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 @@ -27244,7 +27356,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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,& @@ -27257,8 +27369,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! END IF !c!------------------------------------------------------------------- !c! NAPISY KONCOWE - END DO ! j - END DO ! iint + ! END DO ! j + !END DO ! iint END DO ! i !c write (iout,*) "Number of loop steps in EGB:",ind !c energy_dec=.false. @@ -27271,7 +27383,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) @@ -27310,7 +27422,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !c! Coulomb electrostatic interaction Ecl = (332.0d0 * Qij) / Rhead !c! derivative of Ecl is Gcl... - dGCLdR = (-332.0d0 * Qij ) / Rhead_sq + dGCLdR = (-332.0d0 * Qij ) / Rhead_sq*sss_ele_cut dGCLdOM1 = 0.0d0 dGCLdOM2 = 0.0d0 dGCLdOM12 = 0.0d0 @@ -27326,7 +27438,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! -(332.0d0 * Qij *& (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb ) - dGGBdR = dGGBdFGB * dFGBdR + 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" @@ -27347,7 +27459,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !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" @@ -27374,9 +27486,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! * ( 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 + dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut !c! dPOLdR1 = 0.0d0 - dPOLdR2 = dPOLdFGB2 * dFGBdR2 + dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut !c! dPOLdR2 = 0.0d0 dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 !c! dPOLdOM1 = 0.0d0 @@ -27389,7 +27501,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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))) + + (( 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 @@ -27419,7 +27531,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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& @@ -27428,14 +27540,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! - 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)& @@ -27443,7 +27555,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! - 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) & @@ -27451,7 +27563,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! + dGCVdR * erhead(k) & + dPOLdR1 * erhead_tail(k,1) & + dPOLdR2 * erhead_tail(k,2)& - + dGLJdR * erhead(k) + + dGLJdR * erhead(k)+sgrad END DO RETURN @@ -27500,7 +27612,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !c! Coulomb electrostatic interaction Ecl = (332.0d0 * Qij) / Rhead !c! derivative of Ecl is Gcl... - dGCLdR = (-332.0d0 * Qij ) / Rhead_sq + dGCLdR = (-332.0d0 * Qij ) / Rhead_sq*sss_ele_cut+ECL*sss_ele_grad + ECL=ECL*sss_ele_cut dGCLdOM1 = 0.0d0 dGCLdOM2 = 0.0d0 dGCLdOM12 = 0.0d0 @@ -27518,7 +27631,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! -(332.0d0 * Qij *& (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb ) - dGGBdR = dGGBdFGB * dFGBdR + dGGBdR = dGGBdFGB * dFGBdR*sss_ele_cut+Egb*sss_ele_grad + Egb=Egb*sss_ele_grad !c!------------------------------------------------------------------- !c! Fisocav - isotropic cavity creation term !c! or "how much energy it costs to put charged head in water" @@ -27539,7 +27653,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !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& + +FisoCav*sss_ele_grad + FisoCav=FisoCav*sss_ele_cut !c!------------------------------------------------------------------- !c! Epol !c! Polarization energy - charged heads polarize hydrophobic "neck" @@ -27566,13 +27682,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! * ( 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 + dPOLdR1 = dPOLdFGB1 * dFGBdR1!*sss_ele_cut+epol*sss_ele_grad !c! dPOLdR1 = 0.0d0 - dPOLdR2 = dPOLdFGB2 * dFGBdR2 + 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 @@ -27581,7 +27698,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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))) + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut& + +(Elj+epol)*sss_ele_grad + Elj=Elj*sss_ele_cut + epol=epol*sss_ele_cut !c!------------------------------------------------------------------- !c! Return the results !c! These things do the dRdX derivatives, that is @@ -27658,7 +27778,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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" @@ -27699,15 +27819,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k)) gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k) !c! this acts on hydrophobic center of interaction - gvdwx(k,i)= gvdwx(k,i) - gg(k) & + gvdwx(k,i)= gvdwx(k,i) - gg(k)*sss_ele_cut & + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))& - + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv - gvdwx(k,j)= gvdwx(k,j) + gg(k) & + + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss_ele_cut + gvdwx(k,j)= gvdwx(k,j) + gg(k)*sss_ele_cut & + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))& - + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss_ele_cut !c! this acts on Calpha - gvdwc(k,i)=gvdwc(k,i)-gg(k) - gvdwc(k,j)=gvdwc(k,j)+gg(k) + gvdwc(k,i)=gvdwc(k,i)-gg(k)*sss_ele_cut + gvdwc(k,j)=gvdwc(k,j)+gg(k)*sss_ele_cut END DO !c! sc_grad is done, now we will compute eheadtail = 0.0d0 @@ -27814,6 +27934,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2) dbot = 12.0d0 * al4 * pom ** 11.0d0 dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig + !c! dGCVdR = 0.0d0 !c!------------------------------------------------------------------- !c! Polarization energy @@ -27920,6 +28041,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 & @@ -27931,7 +28053,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! - dQUADdR * pom& - tuna(k) & + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))& - + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv + + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss_ele_cut pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) !c! this acts on hydrophobic center of interaction @@ -27946,7 +28068,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! + dQUADdR * pom & + tuna(k) & + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) & - + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss_ele_cut !c! this acts on Calpha gheadtail(k,3,1) = gheadtail(k,3,1) & @@ -27992,16 +28114,22 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -28045,7 +28173,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) & * (2.0d0 - 0.5d0 * ee1) ) & / (2.0d0 * fgb1) - dPOLdR1 = dPOLdFGB1 * dFGBdR1 + dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut +! epol=epol*sss_ele_cut !c! dPOLdR1 = 0.0d0 dPOLdOM1 = 0.0d0 dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 @@ -28062,13 +28191,16 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -28106,7 +28238,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & * (2.0d0 - 0.5d0 * ee2) ) & / (2.0d0 * fgb2) - dPOLdR2 = dPOLdFGB2 * dFGBdR2 + dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut +! epol=epol*sss_ele_cut !c! dPOLdR2 = 0.0d0 dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 !c! dPOLdOM1 = 0.0d0 @@ -28127,14 +28260,18 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -28173,7 +28310,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & * (2.0d0 - 0.5d0 * ee2) ) & / (2.0d0 * fgb2) - dPOLdR2 = dPOLdFGB2 * dFGBdR2 + dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut+epol*sss_ele_grad + epol=epol*sss_ele_cut !c! dPOLdR2 = 0.0d0 dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 !c! dPOLdOM1 = 0.0d0 @@ -28211,7 +28349,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) @@ -28238,8 +28376,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! hawk = w2 * Qi * Qi * (1.0d0 - sqom2) Ecl = sparrow / Rhead**2.0d0 & - hawk / Rhead**4.0d0 - dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 & - + 4.0d0 * hawk / Rhead**5.0d0 + dGCLdR = (- 2.0d0 * sparrow / Rhead**3.0d0 & + + 4.0d0 * hawk / Rhead**5.0d0)*sss_ele_cut !c! dF/dom1 dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0) !c! dF/dom2 @@ -28263,7 +28401,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) & * (2.0d0 - 0.5d0 * ee1) ) & / (2.0d0 * fgb1) - dPOLdR1 = dPOLdFGB1 * dFGBdR1 + dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut !c! dPOLdR1 = 0.0d0 dPOLdOM1 = 0.0d0 dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 @@ -28275,7 +28413,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !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))) + + (( 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) @@ -28292,114 +28430,235 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 END SUBROUTINE eqd - SUBROUTINE edq(Ecl,Elj,Epol) -! IMPLICIT NONE - use comm_momo - use calc_data - double precision facd3, adler,ecl,elj,epol - alphapol2 = alphapol(itypj,itypi) - w1 = wqdip(1,itypi,itypj) - w2 = wqdip(2,itypi,itypj) - pis = sig0head(itypi,itypj) - eps_head = epshead(itypi,itypj) + SUBROUTINE eqd_cat(Ecl,Elj,Epol) + use calc_data + use comm_momo + double precision facd4, federmaus,ecl,elj,epol + alphapol1 = alphapolcat(itypi,itypj) + w1 = wqdipcat(1,itypi,itypj) + w2 = wqdipcat(2,itypi,itypj) + pis = sig0headcat(itypi,itypj) + eps_head = epsheadcat(itypi,itypj) +! eps_head=0.0d0 +! w2=0.0d0 +! alphapol1=0.0d0 !c!------------------------------------------------------------------- -!c! R2 - distance between head of jth side chain and tail of ith sidechain - R2 = 0.0d0 +!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 - R2=R2+(chead(k,2)-ctail(k,1))**2 + R1=R1+(ctail(k,2)-chead(k,1))**2 END DO !c! Pitagoras - R2 = dsqrt(R2) + 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 * Qj * om1 - hawk = w2 * Qj * Qj * (1.0d0 - sqom2) - ECL = sparrow / Rhead**2.0d0 & + sparrow = w1 * Qi * om1 + hawk = w2 * Qi * Qi * (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 + dGCLdR =sss_ele_cut*(-2.0d0 * sparrow / Rhead**3.0d0 & + + 4.0d0 * hawk / Rhead**5.0d0)+sss_ele_grad*ECL + ECL=ECL*sss_ele_cut !c! dF/dom1 - dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0) + dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0) !c! dF/dom2 - dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0) + dGCLdOM2 = 0.0d0 ! + +!(2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0) + !c-------------------------------------------------------------------- !c Polarization energy !c Epol - MomoFac2 = (1.0d0 - chi2 * sqom1) - RR2 = R2 * R2 / MomoFac2 - ee2 = exp(-(RR2 / (4.0d0 * a12sq))) - fgb2 = sqrt(RR2 + a12sq * ee2) - epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 ) - dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) & - / (fgb2 ** 5.0d0) - dFGBdR2 = ( (R2 / MomoFac2) & - * ( 2.0d0 - (0.5d0 * ee2) ) ) & - / (2.0d0 * fgb2) - dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & - * (2.0d0 - 0.5d0 * ee2) ) & - / (2.0d0 * fgb2) - dPOLdR2 = dPOLdFGB2 * dFGBdR2 -!c! dPOLdR2 = 0.0d0 - dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 -!c! dPOLdOM1 = 0.0d0 + 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+epol*sss_ele_grad +!c! dPOLdR1 = 0.0d0 + dPOLdOM1 = 0.0d0 +! dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 dPOLdOM2 = 0.0d0 + epol=epol*sss_ele_cut !c!------------------------------------------------------------------- !c! Elj pom = (pis / Rhead)**6.0d0 Elj = 4.0d0 * eps_head * pom * (pom-1.0d0) !c! derivative of Elj is Glj - dGLJdR = 4.0d0 * eps_head & - * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) & - + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0))) -!c!------------------------------------------------------------------- -!c! Return the results -!c! (see comments in Eqq) + dGLJdR = 4.0d0 * eps_head*sss_ele_cut & + * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) & + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))+Elj*sss_ele_grad + Elj=Elj*sss_ele_cut DO k = 1, 3 erhead(k) = Rhead_distance(k)/Rhead - erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1) END DO + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) - erdxj = scalar( erhead(1), dC_norm(1,j+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)) + gradpepcatx(k,i) = gradpepcatx(k,i) & + - dGCLdR * pom& + - dPOLdR1 * hawk & + - dGLJdR * pom + +! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) +! gradpepcatx(k,j) = gradpepcatx(k,j) & +! + dGCLdR * pom & +! + dPOLdR1 * (erhead_tail(k,1) & +! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) & +! + dGLJdR * pom + + + gradpepcat(k,i) = gradpepcat(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR1 * erhead_tail(k,1) & + - dGLJdR * erhead(k) + + gradpepcat(k,j) = gradpepcat(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) & + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE eqd_cat + + SUBROUTINE edq(Ecl,Elj,Epol) +! IMPLICIT NONE + use comm_momo + use calc_data + + double precision facd3, adler,ecl,elj,epol,sgrad + alphapol2 = alphapol(itypj,itypi) + w1 = wqdip(1,itypi,itypj) + w2 = wqdip(2,itypi,itypj) + pis = sig0head(itypi,itypj) + eps_head = epshead(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) + ECL = sparrow / Rhead**2.0d0 & + - hawk / Rhead**4.0d0 +!c!------------------------------------------------------------------- +!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) +!c! dF/dom1 + dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0) +!c! dF/dom2 + dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0) +!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+nres) ) eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) facd1 = d1 * vbld_inv(i+nres) @@ -28408,30 +28667,30 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -28473,12 +28732,13 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !c!------------------------------------------------------------------- !c! derivative of ecl is Gcl !c! dF/dr part - dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 & - + 4.0d0 * hawk / Rhead**5.0d0 + dGCLdR =( - 2.0d0 * sparrow / Rhead**3.0d0 & + + 4.0d0 * hawk / Rhead**5.0d0)*sss_ele_cut+ECL*sss_ele_grad !c! dF/dom1 dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0) !c! dF/dom2 dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0) + ECL=ECL*sss_ele_cut !c-------------------------------------------------------------------- !c-------------------------------------------------------------------- !c Polarization energy @@ -28496,11 +28756,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & * (2.0d0 - 0.5d0 * ee2) ) & / (2.0d0 * fgb2) - dPOLdR2 = dPOLdFGB2 * dFGBdR2 + dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut+epol*sss_ele_grad !c! dPOLdR2 = 0.0d0 dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 !c! dPOLdOM1 = 0.0d0 dPOLdOM2 = 0.0d0 + epol=epol*sss_ele_cut !c!------------------------------------------------------------------- !c! Elj pom = (pis / Rhead)**6.0d0 @@ -28508,7 +28769,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !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))) + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut+& + Elj*sss_ele_grad + Elj=Elj*sss_ele_cut !c!------------------------------------------------------------------- !c! Return the results @@ -28593,8 +28856,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !c!------------------------------------------------------------------- !c! derivative of ecl is Gcl !c! dF/dr part - dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 & - + 4.0d0 * hawk / Rhead**5.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 !c! dF/dom1 dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0) !c! dF/dom2 @@ -28616,7 +28881,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & * (2.0d0 - 0.5d0 * ee2) ) & / (2.0d0 * fgb2) - dPOLdR2 = dPOLdFGB2 * dFGBdR2 + dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut+epol*sss_ele_grad + epol=epol*sss_ele_grad !c! dPOLdR2 = 0.0d0 dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 !c! dPOLdOM1 = 0.0d0 @@ -28626,9 +28892,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 & + dGLJdR = 4.0d0 * eps_head*sss_ele_cut & * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) & - + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0))) + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))+Elj*sss_ele_grad + Elj=Elj*sss_ele_cut !c!------------------------------------------------------------------- !c! Return the results @@ -28715,7 +28982,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 + 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) & @@ -28743,15 +29010,141 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 + SUBROUTINE edd_cat(ECL) +! IMPLICIT NONE + use comm_momo + use calc_data + + double precision ecl +!c! csig = sigiso(itypi,itypj) + w1 = wqdipcat(1,itypi,itypj) + w2 = wqdipcat(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 simplification + 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+ECL*sss_ele_grad +!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)) + gradpepcatx(k,i) = gradpepcatx(k,i) - dGCLdR * pom +! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) +! gradpepcatx(k,j) = gradpepcatx(k,j) + dGCLdR * pom + + gradpepcat(k,i) = gradpepcat(k,i) - dGCLdR * erhead(k) + gradpepcat(k,j) = gradpepcat(k,j) + dGCLdR * erhead(k) + END DO + RETURN + END SUBROUTINE edd_cat + SUBROUTINE edd_cat_pep(ECL) +! IMPLICIT NONE + use comm_momo + use calc_data + + double precision ecl +!c! csig = sigiso(itypi,itypj) + w1 = wqdipcat(1,itypi,itypj) + w2 = wqdipcat(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+ECL*sss_ele_grad + ECL=ECL*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)) + gradpepcat(k,i) = gradpepcat(k,i) + dGCLdR * pom + gradpepcat(k,i+1) = gradpepcat(k,i+1) - dGCLdR * pom +! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) +! gradpepcatx(k,j) = gradpepcatx(k,j) + dGCLdR * pom + + gradpepcat(k,i) = gradpepcat(k,i) - dGCLdR * erhead(k)*0.5d0 + gradpepcat(k,i+1) = gradpepcat(k,i+1)- dGCLdR * erhead(k)*0.5d0 + gradpepcat(k,j) = gradpepcat(k,j) + dGCLdR * erhead(k) + END DO + RETURN + END SUBROUTINE edd_cat_pep + SUBROUTINE elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol) ! IMPLICIT NONE use comm_momo @@ -29718,7 +30111,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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:: contlistcatpnormi(250*nres),contlistcatpnormj(250*nres) + integer:: contlistcatpnormi(300*nres),contlistcatpnormj(300*nres) integer:: contlistcatscnormi(250*nres),contlistcatscnormj(250*nres) integer:: contlistcatptrani(250*nres),contlistcatptranj(250*nres) integer:: contlistcatsctrani(250*nres),contlistcatsctranj(250*nres) @@ -29764,9 +30157,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -29813,8 +30206,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! if (itype(j,5).le.5) then ilist_catscnorm=ilist_catscnorm+1 ! this can be substituted by cantor and anti-cantor +! write(iout,*) "have contact",i,j,ilist_catscnorm contlistcatscnormi(ilist_catscnorm)=i contlistcatscnormj(ilist_catscnorm)=j +! write(iout,*) "have contact2",i,j,ilist_catscnorm,& +! contlistcatscnormi(ilist_catscnorm),contlistcatscnormj(ilist_catscnorm) + else ilist_catsctran=ilist_catsctran+1 ! this can be substituted by cantor and anti-cantor @@ -29841,16 +30238,17 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ilist_catscnorm,ilist_catpnorm,ilist_catscang do i=1,ilist_catsctran - write (iout,*) i,contlistcatsctrani(i),contlistcatsctranj(i) + 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),contlistcatsctranj(i) + write (iout,*) i,contlistcatscnormi(i),contlistcatscnormj(i) enddo do i=1,ilist_catpnorm - write (iout,*) i,contlistcatpnormi(i),contlistcatsctranj(i) + write (iout,*) i,contlistcatpnormi(i),contlistcatscnormj(i) enddo do i=1,ilist_catscang write (iout,*) i,contlistcatscangi(i),contlistcatscangi(i) @@ -30162,168 +30560,504 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) + -!----------------------------------------------------------------------------- - double precision function boxshift(x,boxsize) - implicit none - double precision x,boxsize - double precision xtemp - xtemp=dmod(x,boxsize) - if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then - boxshift=xtemp-boxsize - else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then - boxshift=xtemp+boxsize - else - boxshift=xtemp - endif - return - end function boxshift -!----------------------------------------------------------------------------- - subroutine to_box(xi,yi,zi) - implicit none -! include 'DIMENSIONS' -! include 'COMMON.CHAIN' - double precision xi,yi,zi - xi=dmod(xi,boxxsize) - if (xi.lt.0.0d0) xi=xi+boxxsize - yi=dmod(yi,boxysize) - if (yi.lt.0.0d0) yi=yi+boxysize - zi=dmod(zi,boxzsize) - if (zi.lt.0.0d0) zi=zi+boxzsize - return - end subroutine to_box -!-------------------------------------------------------------------------- - subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi) - implicit none -! include 'DIMENSIONS' -! include 'COMMON.IOUNITS' -! include 'COMMON.CHAIN' - double precision xi,yi,zi,sslipi,ssgradlipi - double precision fracinbuf -! double precision sscalelip,sscagradlip -#ifdef DEBUG - write (iout,*) "bordlipbot",bordlipbot," bordliptop",bordliptop - write (iout,*) "buflipbot",buflipbot," lipbufthick",lipbufthick - write (iout,*) "xi yi zi",xi,yi,zi -#endif - if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then -! the energy transfer exist - if (zi.lt.buflipbot) then -! what fraction I am in - fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick) -! lipbufthick is thickenes of lipid buffore - sslipi=sscalelip(fracinbuf) - ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) - sslipi=sscalelip(fracinbuf) - ssgradlipi=sscagradlip(fracinbuf)/lipbufthick else - sslipi=1.0d0 - ssgradlipi=0.0 + 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 - else - sslipi=0.0d0 - ssgradlipi=0.0 - 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,*) "sslipi",sslipi," ssgradlipi",ssgradlipi + 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 lipid_layer -!------------------------------------------------------------- - subroutine ecat_prot_transition(ecation_prottran) - integer:: itypi,itypj,ityptrani,ityptranj,k,l,i,j - real(kind=8),dimension(3):: cjtemp,citemp,diff,dsctemp,vecsc,& - diffnorm,boxx,r,dEvan1Cm,dEvan2Cm,dEtotalCm - real(kind=8):: ecation_prottran,dista,sdist,De,ene,x0left,& - alphac,grad,sumvec,simplesum,pom,erdxi,facd1,& - sss_ele_cut,sss_ele_cut_grad,sss2min,sss2mingrad,& - ene1,ene2,grad1,grad2,evan1,evan2,rcal,r4,r7,r0p,& - r06,r012,epscalc,rocal,ract - ecation_prottran=0.0d0 - boxx(1)=boxxsize - boxx(2)=boxysize - boxx(3)=boxzsize - do k=g_listcatsctran_start,g_listcatsctran_end - i=newcontlistcatsctrani(k) - j=newcontlistcatsctranj(k) -! print *,i,j,"in new tran" - do l=1,3 - citemp(l)=c(l,i+nres) - cjtemp(l)=c(l,j) - enddo + end subroutine make_lip_pep_list - itypi=itype(i,1) !as the first is the protein part - itypj=itype(j,5) !as the second part is always cation -! remapping to internal types -! read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),& -! (athetacattran(k,j,i),k=1,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),& -! demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),& -! x0cattrans(j,i) - - if (itypj.eq.6) then - ityptranj=1 !as now only Zn2+ is this needs to be modified for other ions - endif - if (itypi.eq.16) then - ityptrani=1 - elseif (itypi.eq.1) then - ityptrani=2 - elseif (itypi.eq.15) then - ityptrani=3 - elseif (itypi.eq.17) then - ityptrani=4 - elseif (itypi.eq.2) then - ityptrani=5 - else - ityptrani=6 - endif - if (ityptrani.gt.ntrantyp(ityptranj)) then -! do l=1,3 -! write(iout,*),gradcattranc(l,j),gradcattranx(l,i) -! enddo -!volume excluded - call to_box(cjtemp(1),cjtemp(2),cjtemp(3)) - call to_box(citemp(1),citemp(2),citemp(3)) - rcal=0.0d0 - do l=1,3 - r(l)=boxshift(cjtemp(l)-citemp(l),boxx(l)) - rcal=rcal+r(l)*r(l) - enddo - ract=sqrt(rcal) - if (ract.gt.r_cut_ele) cycle - sss_ele_cut=sscale_ele(ract) - sss_ele_cut_grad=sscagrad_ele(ract) - rocal=1.5 - epscalc=0.2 - r0p=0.5*(rocal+sig0(itype(i,1))) - r06 = r0p**6 - r012 = r06*r06 - Evan1=epscalc*(r012/rcal**6) - Evan2=epscalc*2*(r06/rcal**3) - r4 = rcal**4 - r7 = rcal**7 - do l=1,3 - dEvan1Cm(l) = 12*r(l)*epscalc*r012/r7 - dEvan2Cm(l) = 12*r(l)*epscalc*r06/r4 - enddo - do l=1,3 - dEtotalCm(l)=(dEvan1Cm(l)+dEvan2Cm(l))*sss_ele_cut-& - (Evan1+Evan2)*sss_ele_cut_grad*r(l)/ract - enddo - ecation_prottran = ecation_prottran+& - (Evan1+Evan2)*sss_ele_cut - do l=1,3 - gradcattranx(l,i)=gradcattranx(l,i)+dEtotalCm(l) - gradcattranc(l,i)=gradcattranc(l,i)+dEtotalCm(l) - gradcattranc(l,j)=gradcattranc(l,j)-dEtotalCm(l) - enddo + 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 + 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,dimension(:),allocatable:: contlistcatpnormi,contlistcatpnormj +! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres) + integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_catscnorm,& + ilist_catsctran,ilist_catpnorm,ilist_catptran,itmp,ilist_catscang,& + ilist_catscangf,ilist_catscangt,k + integer displ(0:nprocs),i_ilist_catscnorm(0:nprocs),ierr,& + i_ilist_catpnorm(0:nprocs),i_ilist_catsctran(0:nprocs),& + i_ilist_catptran(0:nprocs),i_ilist_catscang(0:nprocs),& + i_ilist_catscangf(0:nprocs),i_ilist_catscangt(0:nprocs) +! write(iout,*),"START make_catcat" + ilist_catpnorm=0 + ilist_catscnorm=0 + ilist_catptran=0 + ilist_catsctran=0 + ilist_catscang=0 - ene=0.0d0 - else -! cycle + if (.not.allocated(contlistcatpnormi)) then + allocate(contlistcatpnormi(900*nres)) + allocate(contlistcatpnormj(900*nres)) + endif + r_buff_list=3.0 + 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=icatb_start,icatb_end + xi=c(1,i) + yi=c(2,i) + zi=c(3,i) + 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+nres) + do j=i+1,itmp+nres_molec(5) + 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-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.(10.0+r_buff_list)) then +! Here the list is created +! if (i.eq.2) then +! print *,i,j,dist_init,ilist_catpnorm +! endif + ilist_catpnorm=ilist_catpnorm+1 + +! this can be substituted by cantor and anti-cantor + contlistcatpnormi(ilist_catpnorm)=i + contlistcatpnormj(ilist_catpnorm)=j + endif +! enddo + enddo + enddo +#ifdef DEBUG + write (iout,*) "before MPIREDUCE",ilist_catsctran,ilist_catptran,& + ilist_catscnorm,ilist_catpnorm,ilist_catscang + + do i=1,ilist_catpnorm + write (iout,*) i,contlistcatpnormi(i) + enddo + + +#endif + if (nfgtasks.gt.1)then + + call MPI_Reduce(ilist_catpnorm,g_ilist_catcatnorm,1,& + MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +! write(iout,*) "before bcast",g_ilist_sc + call MPI_Gather(ilist_catpnorm,1,MPI_INTEGER,& + i_ilist_catpnorm,1,MPI_INTEGER,king,FG_COMM,IERR) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_catpnorm(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlistcatpnormi,ilist_catpnorm,MPI_INTEGER,& + newcontlistcatcatnormi,i_ilist_catpnorm,displ,MPI_INTEGER,& + king,FG_COMM,IERR) + call MPI_Gatherv(contlistcatpnormj,ilist_catpnorm,MPI_INTEGER,& + newcontlistcatcatnormj,i_ilist_catpnorm,displ,MPI_INTEGER,& + king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_catcatnorm,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(newcontlistcatcatnormi,g_ilist_catcatnorm,MPI_INT,king,FG_COMM,IERR) + call MPI_Bcast(newcontlistcatcatnormj,g_ilist_catcatnorm,MPI_INT,king,FG_COMM,IERR) + + + else + g_ilist_catcatnorm=ilist_catpnorm + do i=1,ilist_catpnorm + newcontlistcatcatnormi(i)=contlistcatpnormi(i) + newcontlistcatcatnormj(i)=contlistcatpnormj(i) + enddo + endif + call int_bounds(g_ilist_catcatnorm,g_listcatcatnorm_start,g_listcatcatnorm_end) + +#ifdef DEBUG + write (iout,*) "after MPIREDUCE",g_ilist_catcatnorm + + do i=1,g_ilist_catcatnorm + write (iout,*) i,newcontlistcatcatnormi(i),newcontlistcatcatnormj(i) + enddo +#endif +! write(iout,*),"END make_catcat" + return + end subroutine make_cat_cat_list + + +!----------------------------------------------------------------------------- + double precision function boxshift(x,boxsize) + implicit none + double precision x,boxsize + double precision xtemp + xtemp=dmod(x,boxsize) + if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then + boxshift=xtemp-boxsize + else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then + boxshift=xtemp+boxsize + else + boxshift=xtemp + endif + return + end function boxshift +!----------------------------------------------------------------------------- + subroutine to_box(xi,yi,zi) + implicit none +! include 'DIMENSIONS' +! include 'COMMON.CHAIN' + double precision xi,yi,zi + xi=dmod(xi,boxxsize) + if (xi.lt.0.0d0) xi=xi+boxxsize + yi=dmod(yi,boxysize) + if (yi.lt.0.0d0) yi=yi+boxysize + zi=dmod(zi,boxzsize) + if (zi.lt.0.0d0) zi=zi+boxzsize + return + end subroutine to_box +!-------------------------------------------------------------------------- + subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi) + implicit none +! include 'DIMENSIONS' +! include 'COMMON.IOUNITS' +! include 'COMMON.CHAIN' + double precision xi,yi,zi,sslipi,ssgradlipi + double precision fracinbuf +! double precision sscalelip,sscagradlip +#ifdef DEBUG + write (iout,*) "bordlipbot",bordlipbot," bordliptop",bordliptop + write (iout,*) "buflipbot",buflipbot," lipbufthick",lipbufthick + write (iout,*) "xi yi zi",xi,yi,zi +#endif + if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then +! the energy transfer exist + if (zi.lt.buflipbot) then +! what fraction I am in + fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick) +! lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif +#ifdef DEBUG + write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi +#endif + return + end subroutine lipid_layer +!------------------------------------------------------------- + subroutine ecat_prot_transition(ecation_prottran) + integer:: itypi,itypj,ityptrani,ityptranj,k,l,i,j + real(kind=8),dimension(3):: cjtemp,citemp,diff,dsctemp,vecsc,& + diffnorm,boxx,r,dEvan1Cm,dEvan2Cm,dEtotalCm + real(kind=8):: ecation_prottran,dista,sdist,De,ene,x0left,& + alphac,grad,sumvec,simplesum,pom,erdxi,facd1,& + sss_ele_cut,sss_ele_cut_grad,sss2min,sss2mingrad,& + ene1,ene2,grad1,grad2,evan1,evan2,rcal,r4,r7,r0p,& + r06,r012,epscalc,rocal,ract + ecation_prottran=0.0d0 + boxx(1)=boxxsize + boxx(2)=boxysize + boxx(3)=boxzsize + write(iout,*) "start ecattran",g_listcatsctran_start,g_listcatsctran_end + do k=g_listcatsctran_start,g_listcatsctran_end + i=newcontlistcatsctrani(k) + j=newcontlistcatsctranj(k) +! print *,i,j,"in new tran" + do l=1,3 + citemp(l)=c(l,i+nres) + cjtemp(l)=c(l,j) + enddo + + itypi=itype(i,1) !as the first is the protein part + itypj=itype(j,5) !as the second part is always cation +! remapping to internal types +! read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),& +! (athetacattran(k,j,i),k=1,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),& +! demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),& +! x0cattrans(j,i) + + if (itypj.eq.6) then + ityptranj=1 !as now only Zn2+ is this needs to be modified for other ions + endif + if (itypi.eq.16) then + ityptrani=1 + elseif (itypi.eq.1) then + ityptrani=2 + elseif (itypi.eq.15) then + ityptrani=3 + elseif (itypi.eq.17) then + ityptrani=4 + elseif (itypi.eq.2) then + ityptrani=5 + else + ityptrani=6 + endif + + if (ityptrani.gt.ntrantyp(ityptranj)) then +! do l=1,3 +! write(iout,*),gradcattranc(l,j),gradcattranx(l,i) +! enddo +!volume excluded + call to_box(cjtemp(1),cjtemp(2),cjtemp(3)) + call to_box(citemp(1),citemp(2),citemp(3)) + rcal=0.0d0 + do l=1,3 + r(l)=boxshift(cjtemp(l)-citemp(l),boxx(l)) + rcal=rcal+r(l)*r(l) + enddo + ract=sqrt(rcal) + if (ract.gt.r_cut_ele) cycle + sss_ele_cut=sscale_ele(ract) + sss_ele_cut_grad=sscagrad_ele(ract) + rocal=1.5 + epscalc=0.2 + r0p=0.5*(rocal+sig0(itype(i,1))) + r06 = r0p**6 + r012 = r06*r06 + Evan1=epscalc*(r012/rcal**6) + Evan2=epscalc*2*(r06/rcal**3) + r4 = rcal**4 + r7 = rcal**7 + do l=1,3 + dEvan1Cm(l) = 12*r(l)*epscalc*r012/r7 + dEvan2Cm(l) = 12*r(l)*epscalc*r06/r4 + enddo + do l=1,3 + dEtotalCm(l)=(dEvan1Cm(l)+dEvan2Cm(l))*sss_ele_cut-& + (Evan1+Evan2)*sss_ele_cut_grad*r(l)/ract + enddo + ecation_prottran = ecation_prottran+& + (Evan1+Evan2)*sss_ele_cut + do l=1,3 + gradcattranx(l,i)=gradcattranx(l,i)+dEtotalCm(l) + gradcattranc(l,i)=gradcattranc(l,i)+dEtotalCm(l) + gradcattranc(l,j)=gradcattranc(l,j)-dEtotalCm(l) + enddo + + ene=0.0d0 + else +! cycle sumvec=0.0d0 simplesum=0.0d0 do l=1,3 @@ -31832,5 +32566,1994 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! return end function sq -!-------------------------------------------------------------------------- +#ifdef LBFGS + double precision function funcgrad(x,g) + use MD_data, only: totT,usampl + implicit none + double precision energia(0:n_ene) + double precision x(nvar),g(nvar) + integer i + call var_to_geom(nvar,x) + call zerograd + call chainbuild + call etotal(energia(0)) + call sum_gradient + funcgrad=energia(0) + call cart2intgrad(nvar,g) + if (usampl) then + do i=1,nres-3 + gloc(i,icg)=gloc(i,icg)+dugamma(i) + enddo + do i=1,nres-2 + gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i) + enddo + endif + do i=1,nvar + g(i)=g(i)+gloc(i,icg) + enddo + return + end function funcgrad + subroutine cart2intgrad(n,g) + integer n + double precision g(n) + double precision drt(3,3,nres),rdt(3,3,nres),dp(3,3),& + temp(3,3),prordt(3,3,nres),prodrt(3,3,nres) + double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp + double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,& + cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl + double precision fromto(3,3),aux(6) + integer i,ii,j,jjj,k,l,m,indi,ind,ind1 + logical sideonly + sideonly=.false. + g=0.0d0 + if (sideonly) goto 10 + do i=1,nres-2 + rdt(1,1,i)=-rt(1,2,i) + rdt(1,2,i)= rt(1,1,i) + rdt(1,3,i)= 0.0d0 + rdt(2,1,i)=-rt(2,2,i) + rdt(2,2,i)= rt(2,1,i) + rdt(2,3,i)= 0.0d0 + rdt(3,1,i)=-rt(3,2,i) + rdt(3,2,i)= rt(3,1,i) + rdt(3,3,i)= 0.0d0 + enddo + do i=2,nres-2 + drt(1,1,i)= 0.0d0 + drt(1,2,i)= 0.0d0 + drt(1,3,i)= 0.0d0 + drt(2,1,i)= rt(3,1,i) + drt(2,2,i)= rt(3,2,i) + drt(2,3,i)= rt(3,3,i) + drt(3,1,i)=-rt(2,1,i) + drt(3,2,i)=-rt(2,2,i) + drt(3,3,i)=-rt(2,3,i) + enddo + ind1=0 + do i=1,nres-2 + ind1=ind1+1 + if (n.gt.nphi) then + + do j=1,3 + do k=1,2 + dpjk=0.0D0 + do l=1,3 + dpjk=dpjk+prod(j,l,i)*rdt(l,k,i) + enddo + dp(j,k)=dpjk + prordt(j,k,i)=dp(j,k) + enddo + dp(j,3)=0.0D0 + g(nphi+i)=g(nphi+i)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg) + enddo + xx1(1)=-0.5D0*xloc(2,i+1) + xx1(2)= 0.5D0*xloc(1,i+1) + do j=1,3 + xj=0.0D0 + do k=1,2 + xj=xj+r(j,k,i)*xx1(k) + enddo + xx(j)=xj + enddo + do j=1,3 + rj=0.0D0 + do k=1,3 + rj=rj+prod(j,k,i)*xx(k) + enddo + g(nphi+i)=g(nphi+i)+rj*gradx(j,i+1,icg) + enddo + if (i.lt.nres-2) then + do j=1,3 + dxoiij=0.0D0 + do k=1,3 + dxoiij=dxoiij+dp(j,k)*xrot(k,i+2) + enddo + g(nphi+i)=g(nphi+i)+dxoiij*gradx(j,i+2,icg) + enddo + endif + + endif + + + if (i.gt.1) then + do j=1,3 + do k=1,3 + dpjk=0.0 + do l=2,3 + dpjk=dpjk+prod(j,l,i)*drt(l,k,i) + enddo + dp(j,k)=dpjk + prodrt(j,k,i)=dp(j,k) + enddo + g(i-1)=g(i-1)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg) + enddo + endif + xx(1)= 0.0D0 + xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i) + xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i) + if (i.gt.1) then + do j=1,3 + rj=0.0D0 + do k=2,3 + rj=rj+prod(j,k,i)*xx(k) + enddo + g(i-1)=g(i-1)-rj*gradx(j,i+1,icg) + enddo + endif + if (i.gt.1) then + do j=1,3 + dxoiij=0.0D0 + do k=1,3 + dxoiij=dxoiij+dp(j,k)*xrot(k,i+2) + enddo + g(i-1)=g(i-1)+dxoiij*gradx(j,i+2,icg) + enddo + endif + do j=i+1,nres-2 + ind1=ind1+1 + call build_fromto(i+1,j+1,fromto) + do k=1,3 + do l=1,3 + tempkl=0.0D0 + do m=1,2 + tempkl=tempkl+prordt(k,m,i)*fromto(m,l) + enddo + temp(k,l)=tempkl + enddo + enddo + if (n.gt.nphi) then + do k=1,3 + g(nphi+i)=g(nphi+i)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg) + enddo + do k=1,3 + dxoijk=0.0D0 + do l=1,3 + dxoijk=dxoijk+temp(k,l)*xrot(l,j+2) + enddo + g(nphi+i)=g(nphi+i)+dxoijk*gradx(k,j+2,icg) + enddo + endif + do k=1,3 + do l=1,3 + tempkl=0.0D0 + do m=1,3 + tempkl=tempkl+prodrt(k,m,i)*fromto(m,l) + enddo + temp(k,l)=tempkl + enddo + enddo + if (i.gt.1) then + do k=1,3 + g(i-1)=g(i-1)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg) + enddo + do k=1,3 + dxoijk=0.0D0 + do l=1,3 + dxoijk=dxoijk+temp(k,l)*xrot(l,j+2) + enddo + g(i-1)=g(i-1)+dxoijk*gradx(k,j+2,icg) + enddo + endif + enddo + enddo + + if (nvar.le.nphi+ntheta) return + + 10 continue + do i=2,nres-1 + if (iabs(itype(i,1)).eq.10 .or. itype(i,1).eq.ntyp1& !) cycle + .or. mask_side(i).eq.0 ) cycle + ii=ialph(i,1) + dsci=vbld(i+nres) +#ifdef OSF + alphi=alph(i) + omegi=omeg(i) + if(alphi.ne.alphi) alphi=100.0 + if(omegi.ne.omegi) omegi=-100.0 +#else + alphi=alph(i) + omegi=omeg(i) +#endif + cosalphi=dcos(alphi) + sinalphi=dsin(alphi) + cosomegi=dcos(omegi) + sinomegi=dsin(omegi) + temp(1,1)=-dsci*sinalphi + temp(2,1)= dsci*cosalphi*cosomegi + temp(3,1)=-dsci*cosalphi*sinomegi + temp(1,2)=0.0D0 + temp(2,2)=-dsci*sinalphi*sinomegi + temp(3,2)=-dsci*sinalphi*cosomegi + theta2=pi-0.5D0*theta(i+1) + cost2=dcos(theta2) + sint2=dsin(theta2) + jjj=0 + do j=1,2 + xp=temp(1,j) + yp=temp(2,j) + xxp= xp*cost2+yp*sint2 + yyp=-xp*sint2+yp*cost2 + zzp=temp(3,j) + xx(1)=xxp + xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1) + xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1) + do k=1,3 + dj=0.0D0 + do l=1,3 + dj=dj+prod(k,l,i-1)*xx(l) + enddo + aux(jjj+k)=dj + enddo + jjj=jjj+3 + enddo + do k=1,3 + g(ii)=g(ii)+aux(k)*gradx(k,i,icg) + g(ii+nside)=g(ii+nside)+aux(k+3)*gradx(k,i,icg) + enddo + enddo + return + end subroutine cart2intgrad + + +#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