X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=7d21f897cd2015227fd479ca133f45ecebc18609;hb=540f877fc0c1eaf1389f53eef08fe02d352e25a2;hp=d77030768df48ac79067415e8de7abc3b3a68ee1;hpb=4baa9e481f53d4c89b076f3aece756fc47282649;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index d770307..7d21f89 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -30,7 +30,7 @@ ! Maximum number of SC local term fitting function coefficiants integer,parameter :: maxsccoef=65 ! Maximum number of local shielding effectors - integer,parameter :: maxcontsshi=50 +! integer,parameter :: maxcontsshi=50 !----------------------------------------------------------------------------- ! commom.calc common/calc/ !----------------------------------------------------------------------------- @@ -253,12 +253,12 @@ #ifdef MPI real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw ! shielding effect varibles for MPI -! real(kind=8) fac_shieldbuf(maxres), -! & grad_shield_locbuf(3,maxcontsshi,-1:maxres), -! & grad_shield_sidebuf(3,maxcontsshi,-1:maxres), -! & grad_shieldbuf(3,-1:maxres) -! integer ishield_listbuf(maxres), -! &shield_listbuf(maxcontsshi,maxres) + real(kind=8) fac_shieldbuf(nres), & + grad_shield_locbuf(3,maxcontsshi,-1:nres), & + grad_shield_sidebuf(3,maxcontsshi,-1:nres), & + grad_shieldbuf(3,-1:nres) + integer ishield_listbuf(nres), & + shield_listbuf(maxcontsshi,nres),k,j,i ! print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, ! & " nfgtasks",nfgtasks @@ -393,7 +393,12 @@ ! Gay-Berne potential (shifted LJ, angular dependence). ! 104 call egb(evdw) case (4) +! print *,"MOMO",scelemode + if (scelemode.eq.0) then call egb(evdw) + else + call emomo(evdw) + endif ! goto 107 ! Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). ! 105 call egbv(evdw) @@ -419,6 +424,56 @@ if (shield_mode.eq.2) then call set_shield_fac2 endif + if (nfgtasks.gt.1) then + call MPI_Allgatherv(fac_shield(ivec_start), & + ivec_count(fg_rank1), & + MPI_DOUBLE_PRECISION,fac_shieldbuf(1),ivec_count(0), & + ivec_displ(0), & + MPI_DOUBLE_PRECISION,FG_COMM,IERROR) + call MPI_Allgatherv(shield_list(1,ivec_start), & + ivec_count(fg_rank1), & + MPI_I50,shield_listbuf(1,1),ivec_count(0), & + ivec_displ(0), & + MPI_I50,FG_COMM,IERROR) + call MPI_Allgatherv(ishield_list(ivec_start), & + ivec_count(fg_rank1), & + MPI_INTEGER,ishield_listbuf(1),ivec_count(0), & + ivec_displ(0), & + MPI_INTEGER,FG_COMM,IERROR) + call MPI_Allgatherv(grad_shield(1,ivec_start), & + ivec_count(fg_rank1), & + MPI_UYZ,grad_shieldbuf(1,1),ivec_count(0), & + ivec_displ(0), & + MPI_UYZ,FG_COMM,IERROR) + call MPI_Allgatherv(grad_shield_side(1,1,ivec_start), & + ivec_count(fg_rank1), & + MPI_SHI,grad_shield_sidebuf(1,1,1),ivec_count(0), & + ivec_displ(0), & + MPI_SHI,FG_COMM,IERROR) + call MPI_Allgatherv(grad_shield_loc(1,1,ivec_start), & + ivec_count(fg_rank1), & + MPI_SHI,grad_shield_locbuf(1,1,1),ivec_count(0), & + ivec_displ(0), & + MPI_SHI,FG_COMM,IERROR) + do i=1,nres + fac_shield(i)=fac_shieldbuf(i) + ishield_list(i)=ishield_listbuf(i) + do j=1,3 + grad_shield(j,i)=grad_shieldbuf(j,i) + enddo !j + do j=1,ishield_list(i) + shield_list(j,i)=shield_listbuf(j,i) + do k=1,3 + grad_shield_loc(k,j,i)=grad_shield_locbuf(k,j,i) + grad_shield_side(k,j,i)=grad_shield_sidebuf(k,j,i) + enddo !k + enddo !j + enddo !i + endif + + + + ! print *,"AFTER EGB",ipot,evdw !mc !mc Sep-06: egb takes care of dynamic ss bonds too @@ -432,6 +487,10 @@ #ifdef TIMING time_vec=time_vec+MPI_Wtime()-time01 #endif + + + + ! print *,"Processor",myrank," left VEC_AND_DERIV" if (ipot.lt.6) then #ifdef SPLITELE @@ -636,6 +695,11 @@ call epep_sc_base(epepbase) call eprot_sc_phosphate(escpho) call eprot_pep_phosphate(epeppho) + else + epepbase=0.0 + escbase=0.0 + escpho=0.0 + epeppho=0.0 endif ! call ecatcat(ecationcation) ! print *,"after ebend", ebe_nucl @@ -3149,6 +3213,7 @@ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 & .or. itype(i+3,1).eq.ntyp1 & .or. itype(i+4,1).eq.ntyp1) cycle +! print *,"before2",i,i+3, gshieldc_t4(2,i+3),gshieldc_t4(2,i) dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3191,6 +3256,7 @@ call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) & call eturn4(i,eello_turn4) +! print *,"before",i,i+3, gshieldc_t4(2,i+3),gshieldc_t4(2,i) num_cont_hb(i)=num_conti enddo ! i ! @@ -3953,8 +4019,12 @@ ! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij ! eel_loc_ij=0.0 - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & - 'eelloc',i,j,eel_loc_ij +! if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & +! 'eelloc',i,j,eel_loc_ij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,8f8.3)') & + 'eelloc',i,j,eel_loc_ij,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) +! print *,"EELLOC",i,gel_loc_loc(i-1) + ! if (energy_dec) write (iout,*) "a22",a22," a23",a23," a32",a32," a33",a33 ! if (energy_dec) write (iout,*) "muij",muij ! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) @@ -4050,7 +4120,7 @@ +aggj1(l,4)*muij(4))& *sss_ele_cut & *fac_shield(i)*fac_shield(j) & - *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) !+eel_loc_ij*sss_ele_grad*rmij*xtemp(l) enddo @@ -4531,8 +4601,10 @@ integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,& rlocshield - + j=i+3 +! if (j.ne.20) return +! print *,i,j,gshieldc_t4(2,j),gshieldc_t4(2,j+1) !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! ! Fourth-order contributions @@ -4611,6 +4683,7 @@ iresshield=shield_list(ilist,i) do k=1,3 rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i) +! print *,"rlocshield",rlocshield,grad_shield_side(k,ilist,i),iresshield gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i) @@ -4621,16 +4694,17 @@ do ilist=1,ishield_list(j) iresshield=shield_list(ilist,j) do k=1,3 +! print *,"rlocshieldj",j,rlocshield,grad_shield_side(k,ilist,j),iresshield rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j) gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ & rlocshield & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j) gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) & +rlocshield +! print *,"after", gshieldc_t4(k,iresshield-1),iresshield-1,gshieldc_t4(k,iresshield) enddo enddo - do k=1,3 gshieldc_t4(k,i)=gshieldc_t4(k,i)+ & grad_shield(k,i)*eello_t4/fac_shield(i) @@ -4640,6 +4714,7 @@ grad_shield(k,i)*eello_t4/fac_shield(i) gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+ & grad_shield(k,j)*eello_t4/fac_shield(j) +! print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1) enddo endif @@ -4758,10 +4833,11 @@ call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) +! if (j.lt.nres-1) then gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) - +! endif a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) @@ -4777,10 +4853,17 @@ call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ! write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 +! if (j.lt.nres-1) then +! print *,"juest before",j1, gcorr4_turn(l,j1) gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) - +! if (shield_mode.gt.0) then +! print *,"juest after",j1, gcorr4_turn(l,j1),gshieldc_t4(k,j1),gshieldc_loc_t4(k,j1),gel_loc_turn4(i+2) +! else +! print *,"juest after",j1, gcorr4_turn(l,j1),gel_loc_turn4(i+2) +! endif +! endif enddo gshieldc_t4(3,i)=gshieldc_t4(3,i)+ & ssgradlipi*eello_t4/4.0d0*lipscale @@ -10553,7 +10636,7 @@ wscpho*gvdwc_scpho(j,i)+ & wpeppho*gvdwc_peppho(j,i) - + @@ -10751,7 +10834,11 @@ +0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)& +wvdwpsb*gvdwpsb1(j,i))& +wbond_nucl*gradb_nucl(j,i)+wsbloc*gsbloc(j,i) - +! if (i.eq.21) then +! print *,"in sum",gradc(j,i,icg),wturn4*gcorr4_turn(j,i),& +! wturn4*gshieldc_t4(j,i), & +! wturn4*gshieldc_loc_t4(j,i) +! endif ! if ((i.le.2).and.(i.ge.1)) ! print *,gradc(j,i,icg),& ! gradbufc(j,i),welec*gelc(j,i), & @@ -10857,7 +10944,8 @@ ! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i) enddo - enddo + enddo +!#define DEBUG #ifdef DEBUG write (iout,*) "gloc before adding corr" do i=1,4*nres @@ -10879,6 +10967,7 @@ write (iout,*) i,gloc(i,icg) enddo #endif +!#undef DEBUG #ifdef MPI if (nfgtasks.gt.1) then do j=1,3 @@ -11040,7 +11129,7 @@ endif endif endif -!el#define DEBUG +!#define DEBUG #ifdef DEBUG write (iout,*) "gradc gradx gloc" do i=1,nres @@ -11048,7 +11137,7 @@ i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg) enddo #endif -!el#undef DEBUG +!#undef DEBUG #ifdef TIMING time_sumgradient=time_sumgradient+MPI_Wtime()-time01 #endif @@ -11065,10 +11154,14 @@ ! include 'COMMON.IOUNITS' real(kind=8), dimension(3) :: dcosom1,dcosom2 ! print *,"wchodze" - eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 - eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 + 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 + -2.0D0*alf12*eps3der+sigder*sigsq_om12& + +dCAVdOM12+ dGCLdOM12 ! diagnostics only ! eom1=0.0d0 ! eom2=0.0d0 @@ -11726,6 +11819,7 @@ c(j,k)=c(j,k)+aincr c(j,k+nres)=c(j,k+nres)+aincr enddo + call zerograd call etotal(energia1) etot1=energia1(0) ggg(j)=(etot1-etot)/aincr @@ -11738,6 +11832,7 @@ do j=1,3 c(j,i+nres)=c(j,i+nres)+aincr dc(j,i+nres)=dc(j,i+nres)+aincr + call zerograd call etotal(energia1) etot1=energia1(0) ggg(j+3)=(etot1-etot)/aincr @@ -11790,31 +11885,18 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=1.0D-5 + aincr=1.0D-4 write(iout,*) 'Calling CHECK_ECARTINT.' nf=0 icall=0 - write (iout,*) "Before geom_to_var" call geom_to_var(nvar,x) - write (iout,*) "after geom_to_var" write (iout,*) "split_ene ",split_ene call flush(iout) if (.not.split_ene) then - write(iout,*) 'Calling CHECK_ECARTINT if' + call zerograd call etotal(energia) -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' etot=energia(0) - write (iout,*) "etot",etot - call flush(iout) -!el call enerprint(energia) -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) call cartgrad -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' - write (iout,*) "exit cartgrad" - call flush(iout) icall =1 do i=1,nres write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) @@ -11822,7 +11904,6 @@ do j=1,3 grad_s(j,0)=gcart(j,0) enddo -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) @@ -11830,19 +11911,12 @@ enddo enddo else -write(iout,*) 'Calling CHECK_ECARTIN else.' !- split gradient check call zerograd call etotal_long(energia) !el call enerprint(energia) - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) icall =1 - write (iout,*) "longrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& (gxcart(j,i),j=1,3) @@ -11859,14 +11933,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call zerograd call etotal_short(energia) call enerprint(energia) - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) icall =1 - write (iout,*) "shortrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& (gxcart(j,i),j=1,3) @@ -11902,6 +11970,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) write (iout,*) "ij",i,j," etot1",etot1 @@ -11922,6 +11991,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot2=energia1(0) write (iout,*) "ij",i,j," etot2",etot2 @@ -11953,6 +12023,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) else @@ -11967,7 +12038,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then - call etotal(energia1) + call zerograd + call etotal(energia1) etot2=energia1(0) ggg(j+3)=(etot1-etot2)/(2*aincr) else @@ -12049,12 +12121,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call etotal(energia) etot=energia(0) !el call enerprint(energia) - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) icall =1 do i=1,nres write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) @@ -12065,6 +12132,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) +! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i) + ! if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo @@ -12074,14 +12143,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call zerograd call etotal_long(energia) !el call enerprint(energia) - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) icall =1 - write (iout,*) "longrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& (gxcart(j,i),j=1,3) @@ -12092,20 +12155,15 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) +! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo call zerograd call etotal_short(energia) !el call enerprint(energia) - call flush(iout) - write (iout,*) "enter cartgrad" - call flush(iout) call cartgrad - write (iout,*) "exit cartgrad" - call flush(iout) icall =1 - write (iout,*) "shortrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& (gxcart(j,i),j=1,3) @@ -12141,6 +12199,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #endif ! call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) ! call enerprint(energia1) @@ -12158,6 +12217,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call chainbuild_cart ! call int_from_cart1(.false.) if (.not.split_ene) then + call zerograd call etotal(energia1) etot2=energia1(0) ggg(j)=(etot1-etot2)/(2*aincr) @@ -12188,6 +12248,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) ! write (iout,*) if (.not.split_ene) then + call zerograd call etotal(energia1) etot1=energia1(0) else @@ -12210,6 +12271,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! write (iout,*) "dxnormnormsafe",dsqrt( ! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) if (.not.split_ene) then + call zerograd call etotal(energia1) etot2=energia1(0) ggg(j+3)=(etot1-etot2)/(2*aincr) @@ -14660,7 +14722,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) ! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij - +! print *,"EELLOC",i,gel_loc_loc(i-1) if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij ! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) !d @@ -16231,7 +16293,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! This subrouting calculates total Cartesian coordinate gradient. ! The subroutine chainbuild_cart and energy MUST be called beforehand. ! -!el#define DEBUG +!#define DEBUG #ifdef TIMING time00=MPI_Wtime() #endif @@ -16239,6 +16301,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call sum_gradient #ifdef TIMING #endif +!#define DEBUG !el write (iout,*) "After sum_gradient" #ifdef DEBUG !el write (iout,*) "After sum_gradient" @@ -16247,6 +16310,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) enddo #endif +!#undef DEBUG ! If performing constraint dynamics, add the gradients of the constraint energy if(usampl.and.totT.gt.eq_time) then do i=1,nct @@ -16273,6 +16337,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #endif ! call checkintcartgrad ! write(iout,*) 'calling int_to_cart' +!#define DEBUG #ifdef DEBUG write (iout,*) "gcart, gxcart, gloc before int_to_cart" #endif @@ -16304,6 +16369,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' (gxcart(j,i),j=1,3) enddo #endif +!#undef DEBUG #ifdef CARGRAD #ifdef DEBUG write (iout,*) "CARGRAD" @@ -16333,7 +16399,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #ifdef TIMING time_cartgrad=time_cartgrad+MPI_Wtime()-time00 #endif - !el#undef DEBUG +!#undef DEBUG return end subroutine cartgrad !----------------------------------------------------------------------------- @@ -16628,22 +16694,22 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+ & cost1*dc_norm(j,i-2))/ & vbld(i-1) - domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i) + domicron(j,1,1,i)=-1.0/sint1*dcosomicron(j,1,1,i) dcosomicron(j,1,2,i)=-(dc_norm(j,i-2) & +cost1*(dc_norm(j,i-1+nres)))/ & vbld(i-1+nres) - domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i) + domicron(j,1,2,i)=-1.0/sint1*dcosomicron(j,1,2,i) !C Calculate derivative over second omicron Sci-1,Cai-1 Cai !C Looks messy but better than if in loop dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres) & +cost2*dc_norm(j,i-1))/ & vbld(i) - domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i) + domicron(j,2,1,i)=-1.0/sint2*dcosomicron(j,2,1,i) dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) & +cost2*(-dc_norm(j,i-1+nres)))/ & vbld(i-1+nres) ! write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres) - domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i) + domicron(j,2,2,i)=-1.0/sint2*dcosomicron(j,2,2,i) enddo endif enddo @@ -16704,15 +16770,20 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* & dc_norm(j,i-3))/vbld(i-2) - dphi(j,1,i)=-1/sing*dcosphi(j,1,i) + dphi(j,1,i)=-1.0/sing*dcosphi(j,1,i) dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2* & dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4* & dcostheta(j,1,i) - dphi(j,2,i)=-1/sing*dcosphi(j,2,i) + dphi(j,2,i)=-1.0/sing*dcosphi(j,2,i) dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4* & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp* & dc_norm(j,i-1))/vbld(i) - dphi(j,3,i)=-1/sing*dcosphi(j,3,i) + dphi(j,3,i)=-1.0/sing*dcosphi(j,3,i) +!#define DEBUG +#ifdef DEBUG + write(iout,*) "just after",dphi(j,3,i),sing,dcosphi(j,3,i) +#endif +!#undef DEBUG endif enddo endif @@ -17050,6 +17121,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),& MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,& king,FG_COMM,IERROR) +!#define DEBUG #ifdef DEBUG !d write (iout,*) "Gather dphi" !d call flush(iout) @@ -17058,6 +17130,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3) enddo #endif +!#undef DEBUG call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),& MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,& king,FG_COMM,IERROR) @@ -17075,6 +17148,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #endif endif #endif +!#define DEBUG #ifdef DEBUG write (iout,*) "dtheta after gather" do i=1,nres @@ -17093,6 +17167,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3) enddo #endif +!#undef DEBUG return end subroutine intcartderiv !----------------------------------------------------------------------------- @@ -19394,7 +19469,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo !C now sscale fraction sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield -!C print *,buff_shield,"buff" +! print *,buff_shield,"buff",sh_frac_dist !C now sscale if (sh_frac_dist.le.0.0) cycle !C print *,ishield_list(i),i @@ -19428,6 +19503,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' long=long_r_sidechain(itype(k,1)) costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2) sinthet=short/dist_pep_side*costhet +! print *,"SORT",short,long,sinthet,costhet !C now costhet_grad !C costhet=0.6d0 !C sinthet=0.8 @@ -19477,7 +19553,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo !C print *,sinphi,sinthet VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) & - & /VSolvSphere_div + /VSolvSphere_div !C & *wshield !C now the gradient... do j=1,3 @@ -19499,19 +19575,22 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sinphi/sinthet*costhet*costhet_grad(j)& +sinthet/sinphi*cosphi*cosphi_grad_long(j))) & )*wshield - +! print *, 1.0d0/(-dsqrt(1.0d0-sinphi*sinthet)),& +! sinphi/sinthet,& +! +sinthet/sinphi,"HERE" grad_shield_loc(j,ishield_list(i),i)= & scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*& (1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(& sinthet/sinphi*cosphi*cosphi_grad_loc(j)& ))& *wshield +! print *,grad_shield_loc(j,ishield_list(i),i) enddo VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist enddo fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield) -!C write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i) +! write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i) enddo return end subroutine set_shield_fac2 @@ -23980,14 +24059,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dscj_inv = vbld_inv(j+1)/2.0 ! Gay-berne var's sig0ij = sigma_peppho - chi1=0.0d0 - chi2=0.0d0 +! chi1=0.0d0 +! chi2=0.0d0 chi12 = chi1 * chi2 - chip1=0.0d0 - chip2=0.0d0 +! chip1=0.0d0 +! chip2=0.0d0 chip12 = chip1 * chip2 - chis1 = 0.0d0 - chis2 = 0.0d0 +! chis1 = 0.0d0 +! chis2 = 0.0d0 chis12 = chis1 * chis2 sig1 = sigmap1_peppho sig2 = sigmap2_peppho @@ -24102,4 +24181,1564 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo enddo end subroutine eprot_pep_phosphate +!!!!!!!!!!!!!!!!------------------------------------------------------------- + subroutine emomo(evdw) + use calc_data + use comm_momo +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.GEO' +! include 'COMMON.VAR' +! include 'COMMON.LOCAL' +! include 'COMMON.CHAIN' +! include 'COMMON.DERIV' +! include 'COMMON.NAMES' +! include 'COMMON.INTERACT' +! include 'COMMON.IOUNITS' +! include 'COMMON.CALC' +! include 'COMMON.CONTROL' +! include 'COMMON.SBRIDGE' + logical :: lprn +!el local variables + integer :: iint,itypi1,subchap,isel + real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi + real(kind=8) :: evdw + 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 + 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,& + 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) + eps_out=80.0d0 + sss_ele_cut=1.0d0 +! print *,"EVDW KURW",evdw,nres + do i=iatsc_s,iatsc_e +! print *,"I am in EVDW",i + itypi=iabs(itype(i,1)) +! if (i.ne.47) cycle + if (itypi.eq.ntyp1) cycle + itypi1=iabs(itype(i+1,1)) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + xi=dmod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=dmod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=dmod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize + + if ((zi.gt.bordlipbot) & + .and.(zi.lt.bordliptop)) then +!C the energy transfer exist + if (zi.lt.buflipbot) then +!C what fraction I am in + fracinbuf=1.0d0- & + ((zi-bordlipbot)/lipbufthick) +!C 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 +! print *, sslipi,ssgradlipi + dxi=dc_norm(1,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) +! dsci_inv=dsc_inv(itypi) + dsci_inv=vbld_inv(i+nres) +! write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres) +! write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi +! +! Calculate SC interaction energy. +! + do iint=1,nint_gr(i) + do j=istart(i,iint),iend(i,iint) +! 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) + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & + 'evdw',i,j,evdwij,' ss' +! if (energy_dec) write (iout,*) & +! 'evdw',i,j,evdwij,' ss' + do k=j+1,iend(i,iint) +!C search over all next residues + if (dyn_ss_mask(k)) then +!C check if they are cysteins +!C write(iout,*) 'k=',k + +!c write(iout,*) "PRZED TRI", evdwij +! evdwij_przed_tri=evdwij + call triple_ssbond_ene(i,j,k,evdwij) +!c if(evdwij_przed_tri.ne.evdwij) then +!c write (iout,*) "TRI:", evdwij, evdwij_przed_tri +!c endif + +!c write(iout,*) "PO TRI", evdwij +!C call the energy function that removes the artifical triple disulfide +!C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & + 'evdw',i,j,evdwij,'tss' + endif!dyn_ss_mask(k) + enddo! k + ELSE +!el ind=ind+1 + itypj=iabs(itype(j,1)) + if (itypj.eq.ntyp1) cycle + CALL elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol) + +! if (j.ne.78) cycle +! dscj_inv=dsc_inv(itypj) + dscj_inv=vbld_inv(j+nres) + xj=c(1,j+nres) + yj=c(2,j+nres) + zj=c(3,j+nres) + xj=dmod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=dmod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=dmod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif + dxj = dc_norm( 1, nres+j ) + dyj = dc_norm( 2, nres+j ) + dzj = dc_norm( 3, nres+j ) +! print *,i,j,itypi,itypj +! d1i=0.0d0 +! d1j=0.0d0 +! BetaT = 1.0d0 / (298.0d0 * Rb) +! Gay-berne var's +!1! sig0ij = sigma_scsc( itypi,itypj ) +! chi1=0.0d0 +! chi2=0.0d0 +! chip1=0.0d0 +! chip2=0.0d0 +! not used by momo potential, but needed by sc_angular which is shared +! by all energy_potential subroutines + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 + a12sq = rborn(itypi,itypj) * rborn(itypj,itypi) +! a12sq = a12sq * a12sq +! charge of amino acid itypi is... + chis1 = chis(itypi,itypj) + chis2 = chis(itypj,itypi) + chis12 = chis1 * chis2 + sig1 = sigmap1(itypi,itypj) + sig2 = sigmap2(itypi,itypj) +! write (*,*) "sig1 = ", sig1 +! chis1=0.0 +! chis2=0.0 +! chis12 = chis1 * chis2 +! sig1=0.0 +! sig2=0.0 +! write (*,*) "sig2 = ", sig2 +! alpha factors from Fcav/Gcav + b1cav = alphasur(1,itypi,itypj) +! b1cav=0.0d0 + b2cav = alphasur(2,itypi,itypj) + b3cav = alphasur(3,itypi,itypj) + b4cav = alphasur(4,itypi,itypj) +! used to determine whether we want to do quadrupole calculations + eps_in = epsintab(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 +! dtail(1,itypi,itypj)=0.0 +! dtail(2,itypi,itypj)=0.0 + + DO k = 1, 3 + ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i) + ctail(k,2)=c(k,j+nres)-dtail(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))) + +! write (*,*) "eps_inout_fac = ", eps_inout_fac +!------------------------------------------------------------------- +! tail location and distance calculations + d1 = dhead(1, 1, itypi, itypj) + d2 = dhead(2, 1, itypi, itypj) + + DO k = 1,3 +! location of polar head is computed by taking hydrophobic centre +! and moving by a d1 * dc_norm vector +! see unres publications for very informative images + chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres) + chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres) +! distance +! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) +! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) + Rhead_distance(k) = chead(k,2) - chead(k,1) + END DO +! pitagoras (root of sum of squares) + Rhead = dsqrt( & + (Rhead_distance(1)*Rhead_distance(1)) & + + (Rhead_distance(2)*Rhead_distance(2)) & + + (Rhead_distance(3)*Rhead_distance(3))) +!------------------------------------------------------------------- +! zero everything that should be zero'ed + evdwij = 0.0d0 + ECL = 0.0d0 + Elj = 0.0d0 + Equad = 0.0d0 + Epol = 0.0d0 + Fcav=0.0d0 + eheadtail = 0.0d0 + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 + Fcav = 0.0d0 + dFdR = 0.0d0 + dCAVdOM1 = 0.0d0 + dCAVdOM2 = 0.0d0 + dCAVdOM12 = 0.0d0 + dscj_inv = vbld_inv(j+nres) +! print *,i,j,dscj_inv,dsci_inv +! rij holds 1/(distance of Calpha atoms) + rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj) + rij = dsqrt(rrij) +!---------------------------- + CALL sc_angular +! this should be in elgrad_init but om's are calculated by sc_angular +! which in turn is used by older potentials +! om = omega, sqom = om^2 + sqom1 = om1 * om1 + sqom2 = om2 * om2 + sqom12 = om12 * om12 + +! now we calculate EGB - Gey-Berne +! It will be summed up in evdwij and saved in evdw + sigsq = 1.0D0 / sigsq + sig = sig0ij * dsqrt(sigsq) +! rij_shift = 1.0D0 / rij - sig + sig0ij + rij_shift = Rtail - sig + sig0ij + IF (rij_shift.le.0.0D0) THEN + evdw = 1.0D20 + RETURN + END IF + sigder = -sig * sigsq + rij_shift = 1.0D0 / rij_shift + fac = rij_shift**expon + c1 = fac * fac * aa_aq(itypi,itypj) +! print *,"ADAM",aa_aq(itypi,itypj) + +! c1 = 0.0d0 + c2 = fac * bb_aq(itypi,itypj) +! c2 = 0.0d0 + evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 ) + eps2der = eps3rt * evdwij + eps3der = eps2rt * evdwij +! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij + evdwij = eps2rt * eps3rt * evdwij +!#ifdef TSCSC +! IF (bb_aq(itypi,itypj).gt.0) THEN +! evdw_p = evdw_p + evdwij +! ELSE +! evdw_m = evdw_m + evdwij +! END IF +!#else + evdw = evdw & + + evdwij +!#endif + + c1 = c1 * eps1 * eps2rt**2 * eps3rt**2 + fac = -expon * (c1 + evdwij) * rij_shift + sigder = fac * sigder +! fac = rij * fac +! Calculate distance derivative + gg(1) = fac + gg(2) = fac + gg(3) = fac +! if (b2.gt.0.0) then + fac = chis1 * sqom1 + chis2 * sqom2 & + - 2.0d0 * chis12 * om1 * om2 * om12 +! we will use pom later in Gcav, so dont mess with it! + pom = 1.0d0 - chis1 * chis2 * sqom12 + Lambf = (1.0d0 - (fac / pom)) +! print *,"fac,pom",fac,pom,Lambf + Lambf = dsqrt(Lambf) + sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0) +! print *,"sig1,sig2",sig1,sig2,itypi,itypj +! write (*,*) "sparrow = ", sparrow + Chif = Rtail * sparrow +! print *,"rij,sparrow",rij , 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 +! print *,top,bot,"bot,top",ChiLambf,Chif + 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) +! dFdL = 0.0d0 + 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+nres) ) + facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + DO k = 1, 3 +!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i) +!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) +!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) +!c! & + ( dFdR * pom ) + + gvdwc(k,i) = gvdwc(k,i) & + - (( dFdR + gg(k) ) * ertail(k)) +!c! & - ( dFdR * ertail(k)) + + gvdwc(k,j) = gvdwc(k,j) & + + (( dFdR + gg(k) ) * ertail(k)) +!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 + + isel = iabs(Qi) + iabs(Qj) +! isel=0 + IF (isel.eq.0) THEN +!c! No charges - do nothing + eheadtail = 0.0d0 + + ELSE IF (isel.eq.4) THEN +!c! Calculate dipole-dipole interactions + CALL edd(ecl) + eheadtail = ECL +! eheadtail = 0.0d0 + + ELSE IF (isel.eq.1 .and. iabs(Qi).eq.1) THEN +!c! Charge-nonpolar interactions + CALL eqn(epol) + eheadtail = epol +! eheadtail = 0.0d0 + + ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN +!c! Nonpolar-charge interactions + CALL enq(epol) + eheadtail = epol +! eheadtail = 0.0d0 + + ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN +!c! Charge-dipole interactions + CALL eqd(ecl, elj, epol) + eheadtail = ECL + elj + epol +! eheadtail = 0.0d0 + + ELSE IF (isel.eq.3 .and. icharge(itypi).eq.2) THEN +!c! Dipole-charge interactions + CALL edq(ecl, elj, epol) + eheadtail = ECL + elj + epol +! eheadtail = 0.0d0 + + ELSE IF ((isel.eq.2.and. & + iabs(Qi).eq.1).and. & + nstate(itypi,itypj).eq.1) THEN +!c! Same charge-charge interaction ( +/+ or -/- ) + CALL eqq(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 -/+ ) + 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 + + 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 + END IF +!c!------------------------------------------------------------------- +!c! NAPISY KONCOWE + END DO ! j + END DO ! iint + END DO ! i +!c write (iout,*) "Number of loop steps in EGB:",ind +!c energy_dec=.false. +! print *,"EVDW KURW",evdw,nres + + RETURN + END SUBROUTINE emomo +!C------------------------------------------------------------------------------------ + SUBROUTINE eqq(Ecl,Egb,Epol,Fisocav,Elj) + use calc_data + use comm_momo + real (kind=8) :: facd3, facd4, federmaus, adler,& + Ecl,Egb,Epol,Fisocav,Elj,Fgb +! integer :: k +!c! Epol and Gpol analytical parameters + alphapol1 = alphapol(itypi,itypj) + alphapol2 = alphapol(itypj,itypi) +!c! Fisocav and Gisocav analytical parameters + al1 = alphiso(1,itypi,itypj) + al2 = alphiso(2,itypi,itypj) + al3 = alphiso(3,itypi,itypj) + al4 = alphiso(4,itypi,itypj) + csig = (1.0d0 & + / dsqrt(sigiso1(itypi, itypj)**2.0d0 & + + sigiso2(itypi,itypj)**2.0d0)) +!c! + pis = sig0head(itypi,itypj) + eps_head = epshead(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) + Egb = -(332.0d0 * Qij * eps_inout_fac) / 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 * eps_inout_fac) / (Fgb * 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 +!c! dPOLdR1 = 0.0d0 + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 +!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+nres) ) + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) + federmaus = scalar(erhead_tail(1,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) + facd2 = d2 * vbld_inv(j+nres) + facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + +!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+nres))) + + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) + gvdwx(k,i) = gvdwx(k,i) & + - 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 + + 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 + + gvdwc(k,i) = gvdwc(k,i) & + - dGCLdR * erhead(k)& + - dGGBdR * erhead(k)& + - dGCVdR * erhead(k)& + - dPOLdR1 * erhead_tail(k,1)& + - dPOLdR2 * erhead_tail(k,2)& + - dGLJdR * erhead(k) + + gvdwc(k,j) = gvdwc(k,j) & + + dGCLdR * erhead(k) & + + dGGBdR * erhead(k) & + + dGCVdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) & + + dPOLdR2 * erhead_tail(k,2)& + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE eqq +!c!------------------------------------------------------------------- + SUBROUTINE energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) + use comm_momo + use calc_data + + double precision eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad + double precision ener(4) + double precision dcosom1(3),dcosom2(3) +!c! used in Epol derivatives + double precision facd3, facd4 + double precision federmaus, adler + integer istate,ii,jj + real (kind=8) :: Fgb +! print *,"CALLING EQUAD" +!c! Epol and Gpol analytical parameters + alphapol1 = alphapol(itypi,itypj) + alphapol2 = alphapol(itypj,itypi) +!c! Fisocav and Gisocav analytical parameters + al1 = alphiso(1,itypi,itypj) + al2 = alphiso(2,itypi,itypj) + al3 = alphiso(3,itypi,itypj) + al4 = alphiso(4,itypi,itypj) + csig = (1.0d0 / dsqrt(sigiso1(itypi, itypj)**2.0d0& + + sigiso2(itypi,itypj)**2.0d0)) +!c! + w1 = wqdip(1,itypi,itypj) + w2 = wqdip(2,itypi,itypj) + pis = sig0head(itypi,itypj) + eps_head = epshead(itypi,itypj) +!c! First things first: +!c! We need to do sc_grad's job with GB and Fcav + eom1 = eps2der * eps2rt_om1 & + - 2.0D0 * alf1 * eps3der& + + sigder * sigsq_om1& + + dCAVdOM1 + eom2 = eps2der * eps2rt_om2 & + + 2.0D0 * alf2 * eps3der& + + sigder * sigsq_om2& + + dCAVdOM2 + eom12 = evdwij * eps1_om12 & + + eps2der * eps2rt_om12 & + - 2.0D0 * alf12 * eps3der& + + sigder *sigsq_om12& + + dCAVdOM12 +!c! now some magical transformations to project gradient into +!c! three cartesian vectors + DO k = 1, 3 + dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k)) + dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k)) + 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) & + + (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) & + + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))& + + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv +!c! this acts on Calpha + gvdwc(k,i)=gvdwc(k,i)-gg(k) + gvdwc(k,j)=gvdwc(k,j)+gg(k) + END DO +!c! sc_grad is done, now we will compute + eheadtail = 0.0d0 + eom1 = 0.0d0 + eom2 = 0.0d0 + eom12 = 0.0d0 + DO istate = 1, nstate(itypi,itypj) +!c************************************************************* + IF (istate.ne.1) THEN + IF (istate.lt.3) THEN + ii = 1 + ELSE + ii = 2 + END IF + jj = istate/ii + d1 = dhead(1,ii,itypi,itypj) + d2 = dhead(2,jj,itypi,itypj) + DO k = 1,3 + chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres) + chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres) + 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))) + END IF + 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 + 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) + Ecl = (332.0d0 * Qij) / (Rhead * eps_in) +!c! Ecl = 0.0d0 +!c! write (*,*) "Ecl = ", Ecl +!c! derivative of Ecl is Gcl... + dGCLdR = (-332.0d0 * Qij ) / (Rhead_sq * eps_in) +!c! dGCLdR = 0.0d0 + dGCLdOM1 = 0.0d0 + dGCLdOM2 = 0.0d0 + dGCLdOM12 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Generalised Born Solvent Polarization + ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq)) + Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0) + Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb +!c! Egb = 0.0d0 +!c! write (*,*) "a1*a2 = ", a12sq +!c! write (*,*) "Rhead = ", Rhead +!c! write (*,*) "Rhead_sq = ", Rhead_sq +!c! write (*,*) "ee = ", ee +!c! write (*,*) "Fgb = ", Fgb +!c! write (*,*) "fac = ", eps_inout_fac +!c! write (*,*) "Qij = ", Qij +!c! write (*,*) "Egb = ", Egb +!c! Derivative of Egb is Ggb... +!c! dFGBdR is used by Quad's later... + dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb) + dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )& + / ( 2.0d0 * Fgb ) + dGGBdR = dGGBdFGB * dFGBdR +!c! dGGBdR = 0.0d0 +!c!------------------------------------------------------------------- +!c! Fisocav - isotropic cavity creation term + pom = Rhead * csig + top = al1 * (dsqrt(pom) + al2 * pom - al3) + bot = (1.0d0 + al4 * pom**12.0d0) + botsq = bot * bot + FisoCav = top / bot + dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2) + dbot = 12.0d0 * al4 * pom ** 11.0d0 + dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig +!c! dGCVdR = 0.0d0 +!c!------------------------------------------------------------------- +!c! Polarization energy +!c! Epol + 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 +!c! derivative of Epol is Gpol... + 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 +!c! dPOLdR1 = 0.0d0 + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 + pom = (pis / Rhead)**6.0d0 + Elj = 4.0d0 * eps_head * pom * (pom-1.0d0) +!c! Elj = 0.0d0 +!c! derivative of Elj is Glj + dGLJdR = 4.0d0 * eps_head & + * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) & + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0))) +!c! dGLJdR = 0.0d0 +!c!------------------------------------------------------------------- +!c! Equad + IF (Wqd.ne.0.0d0) THEN + Beta1 = 5.0d0 + 3.0d0 * (sqom12 - 1.0d0) & + - 37.5d0 * ( sqom1 + sqom2 ) & + + 157.5d0 * ( sqom1 * sqom2 ) & + - 45.0d0 * om1*om2*om12 + fac = -( Wqd / (2.0d0 * Fgb**5.0d0) ) + Equad = fac * Beta1 +!c! Equad = 0.0d0 +!c! derivative of Equad... + dQUADdR = ((2.5d0 * Wqd * Beta1) / (Fgb**6.0d0)) * dFGBdR +!c! dQUADdR = 0.0d0 + dQUADdOM1 = fac* (-75.0d0*om1 + 315.0d0*om1*sqom2 - 45.0d0*om2*om12) +!c! dQUADdOM1 = 0.0d0 + dQUADdOM2 = fac* (-75.0d0*om2 + 315.0d0*sqom1*om2 - 45.0d0*om1*om12) +!c! dQUADdOM2 = 0.0d0 + dQUADdOM12 = fac * ( 6.0d0*om12 - 45.0d0*om1*om2 ) + ELSE + Beta1 = 0.0d0 + Equad = 0.0d0 + END IF +!c!------------------------------------------------------------------- +!c! Return the results +!c! Angular stuff + eom1 = dPOLdOM1 + dQUADdOM1 + eom2 = dPOLdOM2 + dQUADdOM2 + eom12 = dQUADdOM12 +!c! now some magical transformations to project gradient into +!c! three cartesian vectors + DO k = 1, 3 + dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k)) + dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k)) + tuna(k) = eom1 * dcosom1(k) + eom2 * dcosom2(k) + END DO +!c! Radial stuff + 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+nres) ) + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) + federmaus = scalar(erhead_tail(1,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) + facd2 = d2 * vbld_inv(j+nres) + facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres) + facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + 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+nres)) + + pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres)) +!c! this acts on hydrophobic center of interaction + gheadtail(k,1,1) = gheadtail(k,1,1) & + - 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 & + - 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 + + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) +!c! this acts on hydrophobic center of interaction + gheadtail(k,2,1) = gheadtail(k,2,1) & + + 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 & + + 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 + +!c! this acts on Calpha + gheadtail(k,3,1) = gheadtail(k,3,1) & + - dGCLdR * erhead(k)& + - dGGBdR * erhead(k)& + - dGCVdR * erhead(k)& + - dPOLdR1 * erhead_tail(k,1)& + - dPOLdR2 * erhead_tail(k,2)& + - dGLJdR * erhead(k) & + - dQUADdR * erhead(k)& + - tuna(k) +!c! this acts on Calpha + gheadtail(k,4,1) = gheadtail(k,4,1) & + + dGCLdR * erhead(k) & + + dGGBdR * erhead(k) & + + dGCVdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) & + + dPOLdR2 * erhead_tail(k,2) & + + dGLJdR * erhead(k) & + + dQUADdR * erhead(k)& + + tuna(k) + END DO + ener(istate) = ECL + Egb + Epol + Fisocav + Elj + Equad + eheadtail = eheadtail & + + wstate(istate, itypi, itypj) & + * dexp(-betaT * ener(istate)) +!c! foreach cartesian dimension + DO k = 1, 3 +!c! foreach of two gvdwx and gvdwc + DO l = 1, 4 + gheadtail(k,l,2) = gheadtail(k,l,2) & + + wstate( istate, itypi, itypj ) & + * dexp(-betaT * ener(istate)) & + * gheadtail(k,l,1) + gheadtail(k,l,1) = 0.0d0 + END DO + END DO + END DO +!c! Here ended the gigantic DO istate = 1, 4, which starts +!c! at the beggining of the subroutine + + DO k = 1, 3 + 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) + DO l = 1, 4 + gheadtail(k,l,1) = 0.0d0 + gheadtail(k,l,2) = 0.0d0 + END DO + END DO + eheadtail = (-dlog(eheadtail)) / betaT + dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 + dQUADdOM1 = 0.0d0 + dQUADdOM2 = 0.0d0 + dQUADdOM12 = 0.0d0 + RETURN + END SUBROUTINE energy_quad +!!----------------------------------------------------------- + SUBROUTINE eqn(Epol) + use comm_momo + use calc_data + + double precision facd4, federmaus,epol + alphapol1 = alphapol(itypi,itypj) +!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 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) + dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) & + / (fgb1 ** 5.0d0) + dFGBdR1 = ( (R1 / MomoFac1) & + * ( 2.0d0 - (0.5d0 * ee1) ) ) & + / ( 2.0d0 * fgb1 ) + dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) & + * (2.0d0 - 0.5d0 * ee1) ) & + / (2.0d0 * fgb1) + dPOLdR1 = dPOLdFGB1 * dFGBdR1 +!c! dPOLdR1 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 + DO k = 1, 3 + erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1) + END DO + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) + federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres)) + facd1 = d1 * vbld_inv(i+nres) + facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres) + + DO k = 1, 3 + hawk = (erhead_tail(k,1) + & + facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))) + + gvdwx(k,i) = gvdwx(k,i) & + - dPOLdR1 * hawk + gvdwx(k,j) = gvdwx(k,j) & + + dPOLdR1 * (erhead_tail(k,1) & + -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) + + gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1) + gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1) + + END DO + RETURN + END SUBROUTINE eqn + SUBROUTINE enq(Epol) + use calc_data + use comm_momo + double precision facd3, adler,epol + alphapol2 = alphapol(itypj,itypi) +!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 Polarization energy + MomoFac2 = (1.0d0 - chi2 * sqom1) + RR2 = R2 * R2 / MomoFac2 + ee2 = exp(-(RR2 / (4.0d0 * a12sq))) + fgb2 = sqrt(RR2 + a12sq * ee2) + epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 ) + dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) & + / (fgb2 ** 5.0d0) + dFGBdR2 = ( (R2 / MomoFac2) & + * ( 2.0d0 - (0.5d0 * ee2) ) ) & + / (2.0d0 * fgb2) + dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) & + * (2.0d0 - 0.5d0 * ee2) ) & + / (2.0d0 * fgb2) + dPOLdR2 = dPOLdFGB2 * dFGBdR2 +!c! dPOLdR2 = 0.0d0 + dPOLdOM1 = dPOLdFGB2 * dFGBdOM1 +!c! dPOLdOM1 = 0.0d0 + dPOLdOM2 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Return the results +!c! (See comments in Eqq) + DO k = 1, 3 + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) ) + adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) ) + facd2 = d2 * vbld_inv(j+nres) + facd3 = dtail(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+nres))) + + gvdwx(k,i) = gvdwx(k,i) & + - dPOLdR2 * (erhead_tail(k,2) & + -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) + gvdwx(k,j) = gvdwx(k,j) & + + dPOLdR2 * condor + + gvdwc(k,i) = gvdwc(k,i) & + - dPOLdR2 * erhead_tail(k,2) + gvdwc(k,j) = gvdwc(k,j) & + + dPOLdR2 * erhead_tail(k,2) + + END DO + RETURN + END SUBROUTINE enq + SUBROUTINE eqd(Ecl,Elj,Epol) + use calc_data + use comm_momo + double precision facd4, federmaus,ecl,elj,epol + alphapol1 = alphapol(itypi,itypj) + w1 = wqdip(1,itypi,itypj) + w2 = wqdip(2,itypi,itypj) + pis = sig0head(itypi,itypj) + eps_head = epshead(itypi,itypj) +!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 = - 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 = (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 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) & + * (2.0d0 - 0.5d0 * ee1) ) & + / (2.0d0 * fgb1) + dPOLdR1 = dPOLdFGB1 * dFGBdR1 +!c! dPOLdR1 = 0.0d0 + dPOLdOM1 = 0.0d0 + dPOLdOM2 = dPOLdFGB1 * dFGBdOM2 +!c! 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))) + 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) ) + erdxj = scalar( erhead(1), dC_norm(1,j+nres) ) + bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) ) + federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres)) + facd1 = d1 * vbld_inv(i+nres) + facd2 = d2 * vbld_inv(j+nres) + facd4 = dtail(2,itypi,itypj) * vbld_inv(j+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)) + gvdwx(k,i) = gvdwx(k,i) & + - dGCLdR * pom& + - dPOLdR1 * hawk & + - dGLJdR * pom + + 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 + + + gvdwc(k,i) = gvdwc(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR1 * erhead_tail(k,1) & + - dGLJdR * erhead(k) + + gvdwc(k,j) = gvdwc(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR1 * erhead_tail(k,1) & + + dGLJdR * erhead(k) + + 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) +!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 * 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 +!c! dF/dom1 + dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0) +!c! dF/dom2 + dGCLdOM2 = (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 + dPOLdOM2 = 0.0d0 +!c!------------------------------------------------------------------- +!c! Elj + pom = (pis / Rhead)**6.0d0 + Elj = 4.0d0 * eps_head * pom * (pom-1.0d0) +!c! derivative of Elj is Glj + dGLJdR = 4.0d0 * eps_head & + * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) & + + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0))) +!c!------------------------------------------------------------------- +!c! Return the results +!c! (see comments in Eqq) + DO k = 1, 3 + erhead(k) = Rhead_distance(k)/Rhead + erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2) + END DO + erdxi = scalar( erhead(1), dC_norm(1,i+nres) ) + erdxj = scalar( erhead(1), dC_norm(1,j+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) + facd2 = d2 * vbld_inv(j+nres) + facd3 = dtail(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+nres))) + + 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 + + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) + gvdwx(k,j) = gvdwx(k,j) & + + dGCLdR * pom & + + dPOLdR2 * condor & + + dGLJdR * pom + + + gvdwc(k,i) = gvdwc(k,i) & + - dGCLdR * erhead(k) & + - dPOLdR2 * erhead_tail(k,2) & + - dGLJdR * erhead(k) + + gvdwc(k,j) = gvdwc(k,j) & + + dGCLdR * erhead(k) & + + dPOLdR2 * erhead_tail(k,2) & + + dGLJdR * erhead(k) + + END DO + RETURN + END SUBROUTINE edq + SUBROUTINE edd(ECL) +! IMPLICIT NONE + use comm_momo + use calc_data + + double precision ecl +!c! csig = sigiso(itypi,itypj) + w1 = wqdip(1,itypi,itypj) + w2 = wqdip(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! write (*,*) "w1 = ", w1 +!c! write (*,*) "w2 = ", w2 +!c! write (*,*) "om1 = ", om1 +!c! write (*,*) "om2 = ", om2 +!c! write (*,*) "om12 = ", om12 +!c! write (*,*) "fac = ", fac +!c! write (*,*) "c1 = ", c1 +!c! write (*,*) "c2 = ", c2 +!c! write (*,*) "Ecl = ", Ecl +!c! write (*,*) "c2_1 = ", (w2 / Rhead ** 6.0d0) +!c! write (*,*) "c2_2 = ", +!c! & (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2)) +!c!------------------------------------------------------------------- +!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 + fac * fac - 3.0d0 * (sqom1 + sqom2)) + dGCLdR = c1 - c2 +!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 +!c! dECL/dom12 + c1 = w1 / (Rhead ** 3.0d0) + c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0 + 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)) + gvdwx(k,i) = gvdwx(k,i) - dGCLdR * pom + pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres)) + gvdwx(k,j) = gvdwx(k,j) + dGCLdR * pom + + gvdwc(k,i) = gvdwc(k,i) - dGCLdR * erhead(k) + gvdwc(k,j) = gvdwc(k,j) + dGCLdR * erhead(k) + END DO + RETURN + END SUBROUTINE edd + SUBROUTINE elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol) +! IMPLICIT NONE + 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,1) +!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 = sigma( itypi,itypj ) + chi1 = chi( itypi, itypj ) + chi2 = chi( itypj, itypi ) + chi12 = chi1 * chi2 + chip1 = chipp( itypi, itypj ) + chip2 = chipp( itypj, itypi ) + chip12 = chip1 * chip2 +! chi1=0.0 +! chi2=0.0 +! chi12=0.0 +! chip1=0.0 +! chip2=0.0 +! chip12=0.0 +!c! not used by momo potential, but needed by sc_angular which is shared +!c! by all energy_potential subroutines + alf1 = 0.0d0 + alf2 = 0.0d0 + alf12 = 0.0d0 +!c! location, location, location +! xj = c( 1, nres+j ) - xi +! yj = c( 2, nres+j ) - yi +! zj = c( 3, nres+j ) - zi + dxj = dc_norm( 1, nres+j ) + dyj = dc_norm( 2, nres+j ) + dzj = dc_norm( 3, nres+j ) +!c! distance from center of chain(?) to polar/charged head +!c! write (*,*) "istate = ", 1 +!c! write (*,*) "ii = ", 1 +!c! write (*,*) "jj = ", 1 + d1 = dhead(1, 1, itypi, itypj) + d2 = dhead(2, 1, itypi, itypj) +!c! ai*aj from Fgb + a12sq = rborn(itypi,itypj) * rborn(itypj,itypi) +!c! a12sq = a12sq * a12sq +!c! charge of amino acid itypi is... + Qi = icharge(itypi) + Qj = icharge(itypj) + Qij = Qi * Qj +!c! chis1,2,12 + chis1 = chis(itypi,itypj) + chis2 = chis(itypj,itypi) + chis12 = chis1 * chis2 + sig1 = sigmap1(itypi,itypj) + sig2 = sigmap2(itypi,itypj) +!c! write (*,*) "sig1 = ", sig1 +!c! write (*,*) "sig2 = ", sig2 +!c! alpha factors from Fcav/Gcav + b1cav = alphasur(1,itypi,itypj) +! b1cav=0.0 + b2cav = alphasur(2,itypi,itypj) + b3cav = alphasur(3,itypi,itypj) + b4cav = alphasur(4,itypi,itypj) + wqd = wquad(itypi, itypj) +!c! used by Fgb + eps_in = epsintab(itypi,itypj) + eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out)) +!c! write (*,*) "eps_inout_fac = ", eps_inout_fac +!c!------------------------------------------------------------------- +!c! tail location and distance calculations + Rtail = 0.0d0 + DO k = 1, 3 + ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i) + ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j) + END DO +!c! tail distances will be themselves usefull elswhere +!c1 (in Gcav, for example) + Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 ) + Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 ) + Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 ) + Rtail = dsqrt( & + (Rtail_distance(1)*Rtail_distance(1)) & + + (Rtail_distance(2)*Rtail_distance(2)) & + + (Rtail_distance(3)*Rtail_distance(3))) +!c!------------------------------------------------------------------- +!c! Calculate location and distance between polar heads +!c! distance between heads +!c! for each one of our three dimensional space... + d1 = dhead(1, 1, itypi, itypj) + d2 = dhead(2, 1, itypi, itypj) + + DO k = 1,3 +!c! location of polar head is computed by taking hydrophobic centre +!c! and moving by a d1 * dc_norm vector +!c! see unres publications for very informative images + chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres) + chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres) +!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 end module energy