X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=d77030768df48ac79067415e8de7abc3b3a68ee1;hb=4baa9e481f53d4c89b076f3aece756fc47282649;hp=a5d1496d099ec1018335e8d668e4ce3f685c953d;hpb=aeaa0a995cbe76e90e4701f7df3aa7cdd17c8555;p=unres4.git

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index a5d1496..d770307 100644
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -604,7 +604,9 @@
        etube=0.0d0
       endif
 !--------------------------------------------------------
+!      write (iout,*) "NRES_MOLEC(2),",nres_molec(2)
 !      print *,"before",ees,evdw1,ecorr
+      if (nres_molec(2).gt.0) then
       call ebond_nucl(estr_nucl)
       call ebend_nucl(ebe_nucl)
       call etor_nucl(etors_nucl)
@@ -613,12 +615,28 @@
       call epsb(evdwpsb,eelpsb)
       call esb(esbloc)
       call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
+      else
+       etors_nucl=0.0d0
+       estr_nucl=0.0d0
+       ebe_nucl=0.0d0
+       evdwsb=0.0d0
+       eelsb=0.0d0
+       esbloc=0.0d0
+      endif
+      if (nfgtasks.gt.1) then
+      if (fg_rank.eq.0) then
+      call ecatcat(ecationcation)
+      endif
+      else
       call ecatcat(ecationcation)
+      endif
       call ecat_prot(ecation_prot)
+      if (nres_molec(2).gt.0) then
       call eprot_sc_base(escbase)
       call epep_sc_base(epepbase)
       call eprot_sc_phosphate(escpho)
       call eprot_pep_phosphate(epeppho)
+      endif
 !      call ecatcat(ecationcation)
 !      print *,"after ebend", ebe_nucl
 #ifdef TIMING
@@ -2561,13 +2579,21 @@
         endif
 !        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
+           if (itype(i-2,1).eq.0) then
+          iti=ntortyp+1
+           else
           iti = itortyp(itype(i-2,1))
+           endif
         else
           iti=ntortyp+1
         endif
 !        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         if (i.gt. nnt+1 .and. i.lt.nct+1) then
+           if (itype(i-1,1).eq.0) then
+          iti1=ntortyp+1
+           else
           iti1 = itortyp(itype(i-1,1))
+           endif
         else
           iti1=ntortyp+1
         endif
@@ -2608,7 +2634,9 @@
         enddo
 !        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         if (i.gt. nnt+1 .and. i.lt.nct+1) then
-          if (itype(i-1,1).le.ntyp) then
+          if (itype(i-1,1).eq.0) then
+           iti1=ntortyp+1
+          elseif (itype(i-1,1).le.ntyp) then
             iti1 = itortyp(itype(i-1,1))
           else
             iti1=ntortyp+1
@@ -6759,12 +6787,12 @@
         difi=phii-phi0(i)
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else if (difi.lt.-drange(i)) then
           difi=difi+drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         endif
 !        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
 !     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
@@ -6870,12 +6898,12 @@
         difi=pinorm(phii-phi0(i))
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else if (difi.lt.-drange(i)) then
           difi=difi+drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else
           difi=0.0
         endif
@@ -20779,8 +20807,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             evdwij=evdwij*eps2rt
             evdwsb=evdwsb+evdwij
             if (lprn) then
-            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+            sigm=dabs(aa_nucl(itypi,itypj)/bb_nucl(itypi,itypj))**(1.0D0/6.0D0)
+            epsi=bb_nucl(itypi,itypj)**2/aa_nucl(itypi,itypj)
             write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
              restyp(itypi,2),i,restyp(itypj,2),j, &
              epsi,sigm,chi1,chi2,chip1,chip2, &
@@ -21737,14 +21765,17 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         r012 = r06**2
         k0 = 332.0*(2.0*2.0)/80.0
         itmp=0
+        
         do i=1,4
         itmp=itmp+nres_molec(i)
         enddo
-        do i=itmp+1,itmp+nres_molec(i)-1
+!        write(iout,*) "itmp",itmp
+        do i=itmp+1,itmp+nres_molec(5)-1
        
         xi=c(1,i)
         yi=c(2,i)
         zi=c(3,i)
+         
           xi=mod(xi,boxxsize)
           if (xi.lt.0) xi=xi+boxxsize
           yi=mod(yi,boxysize)
@@ -21763,6 +21794,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           if (yj.lt.0) yj=yj+boxysize
           zj=dmod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
+!          write(iout,*) c(1,i),xi,xj,"xy",boxxsize
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -21820,6 +21852,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
         enddo
 
+!        write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
         ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
        enddo
        enddo
@@ -21869,7 +21902,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         do i=1,4
         itmp=itmp+nres_molec(i)
         enddo
-        do i=1,nres_molec(1)-1  ! loop over all peptide groups needs parralelization
+!        do i=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))
@@ -21976,7 +22010,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        enddo ! j
        enddo ! i
 !------------------------------------------sidechains
-        do i=1,nres_molec(1)
+!        do i=1,nres_molec(1)
+        do i=ibond_start,ibond_end
          if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms
 !         cycle
 !        print *,i,ecation_prot
@@ -22409,7 +22444,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        integer troll
        eps_out=80.0d0
        escbase=0.0d0
-       do i=1,nres_molec(1)
+!       do i=1,nres_molec(1)
+        do i=ibond_start,ibond_end
         if (itype(i,1).eq.ntyp1_molec(1)) cycle
         itypi  = itype(i,1)
         dxi    = dc_norm(1,nres+i)
@@ -22944,7 +22980,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        integer troll
        eps_out=80.0d0
        epepbase=0.0d0
-       do i=1,nres_molec(1)-1
+!       do i=1,nres_molec(1)-1
+        do i=ibond_start,ibond_end
         if (itype(i,1).eq.ntyp1_molec(1).or.itype(i+1,1).eq.ntyp1_molec(1)) cycle
 !C        itypi  = itype(i,1)
         dxi    = dc_norm(1,i)
@@ -23346,7 +23383,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       real(kind=8) :: evdw,sig0ij
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
                     dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
-                    sslipi,sslipj,faclip
+                    sslipi,sslipj,faclip,alpha_sco
       integer :: ii
       real(kind=8) :: fracinbuf
        real (kind=8) :: escpho
@@ -23364,7 +23401,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        integer troll
        eps_out=80.0d0
        escpho=0.0d0
-       do i=1,nres_molec(1)
+!       do i=1,nres_molec(1)
+        do i=ibond_start,ibond_end
         if (itype(i,1).eq.ntyp1_molec(1)) cycle
         itypi  = itype(i,1)
         dxi    = dc_norm(1,nres+i)
@@ -23401,6 +23439,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           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
@@ -23627,12 +23666,15 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !c!      write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
 !c!      write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
 !      alphapol1 = alphapol_scpho(itypi)
-       if (wqq_scpho(itypi).gt.0.0) then
+       if (wqq_scpho(itypi).ne.0.0) then
        Qij=wqq_scpho(itypi)/eps_in
+       alpha_sco=1.d0/alphi_scpho(itypi)
 !       Qij=0.0
-       Ecl = (332.0d0 * Qij) / Rhead
+       Ecl = (332.0d0 * Qij*dexp(-Rhead*alpha_sco)) / Rhead
 !c! derivative of Ecl is Gcl...
-       dGCLdR = (-332.0d0 * Qij ) / Rhead_sq
+       dGCLdR = (-332.0d0 * Qij*dexp(-Rhead*alpha_sco)*  &
+                (Rhead*alpha_sco+1) ) / Rhead_sq
+       if (energy_dec) write(iout,*) "ECL",ECL,Rhead,1.0/rij
        else if (wqdip_scpho(2,itypi).gt.0.0d0) then
        w1        = wqdip_scpho(1,itypi)
        w2        = wqdip_scpho(2,itypi)
@@ -23654,6 +23696,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        Ecl = sparrow / Rhead**2.0d0 &
            - hawk    / Rhead**4.0d0
 !c!-------------------------------------------------------------------
+       if (energy_dec) write(iout,*) "ECLdipdip",ECL,Rhead,&
+           1.0/rij,sparrow
+
 !c! derivative of ecl is Gcl
 !c! dF/dr part
        dGCLdR  = - 2.0d0 * sparrow / Rhead**3.0d0 &
@@ -23755,6 +23800,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 
        END DO
 !       if (i.eq.3) print *,i,j,evdwij,epol,Fcav,ECL
+       if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+        "escpho:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,escpho
        escpho=escpho+evdwij+epol+Fcav+ECL
        call sc_grad_scpho
          enddo
@@ -23864,7 +23911,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        integer troll
        real (kind=8) :: dcosom1(3),dcosom2(3)
        epeppho=0.0d0
-       do i=1,nres_molec(1)
+!       do i=1,nres_molec(1)
+        do i=ibond_start,ibond_end
         if (itype(i,1).eq.ntyp1_molec(1)) cycle
         itypi  = itype(i,1)
         dsci_inv = vbld_inv(i+1)/2.0
@@ -23901,6 +23949,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           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