projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
rozgrzebany shield
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 5402971..5cb9a2a 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -2780,6 +2780,17 @@ c       write(iout,*)  'b1=',b1(1,i-2)
 c       write (iout,*) 'theta=', theta(i-1)
        enddo
 #else
+        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          iti = itortyp(itype(i-2))
+        else
+          iti=ntortyp+1
+        endif
+c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+          iti1 = itortyp(itype(i-1))
+        else
+          iti1=ntortyp+1
+        endif
         b1(1,i-2)=b(3,iti)
         b1(2,i-2)=b(5,iti)
         b2(1,i-2)=b(2,iti)
@@ -2934,6 +2945,7 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         do k=1,2
           mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
         enddo
+C        write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
 c        write (iout,*) 'mu ',mu(:,i-2),i-2
 cd        write (iout,*) 'mu1',mu1(:,i-2)
 cd        write (iout,*) 'mu2',mu2(:,i-2)
@@ -3696,8 +3708,8 @@ C 12/26/95 - for the evaluation of multi-body H-bonding interactions
           if (shield_mode.gt.0) then
 C          fac_shield(i)=0.4
 C          fac_shield(j)=0.6
-          el1=el1*fac_shield(i)*fac_shield(j)
-          el2=el2*fac_shield(i)*fac_shield(j)
+          el1=el1*fac_shield(i)**2*fac_shield(j)**2
+          el2=el2*fac_shield(i)**2*fac_shield(j)**2
           eesij=(el1+el2)
           ees=ees+eesij
           else
@@ -3716,7 +3728,8 @@ cd     &      xmedi,ymedi,zmedi,xj,yj,zj
               write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') 
      &'evdw1',i,j,evdwij
      &,iteli,itelj,aaa,evdw1
-              write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+              write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
+     &fac_shield(i),fac_shield(j)
           endif
 
 C
@@ -3743,9 +3756,10 @@ C          print *,i,j
            iresshield=shield_list(ilist,i)
            do k=1,3
            rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
+     &      *2.0
            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
      &              rlocshield
-     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
             gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
 C           gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
 C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
@@ -3768,9 +3782,10 @@ C              endif
            iresshield=shield_list(ilist,j)
            do k=1,3
            rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+     &     *2.0
            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
      &              rlocshield
-     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
 
 C     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
@@ -3793,13 +3808,13 @@ C              endif
 
           do k=1,3
             gshieldc(k,i)=gshieldc(k,i)+
-     &              grad_shield(k,i)*eesij/fac_shield(i)
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
             gshieldc(k,j)=gshieldc(k,j)+
-     &              grad_shield(k,j)*eesij/fac_shield(j)
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
             gshieldc(k,i-1)=gshieldc(k,i-1)+
-     &              grad_shield(k,i)*eesij/fac_shield(i)
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
             gshieldc(k,j-1)=gshieldc(k,j-1)+
-     &              grad_shield(k,j)*eesij/fac_shield(j)
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
 
            enddo
            endif
@@ -3918,7 +3933,7 @@ cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
 cd   &          (dcosg(k),k=1,3)
           do k=1,3
             ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
-     &      fac_shield(i)*fac_shield(j)
+     &      fac_shield(i)**2*fac_shield(j)**2
           enddo
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
@@ -3939,11 +3954,11 @@ C                     print *,"before22", gelc_long(1,i), gelc_long(1,j)
             gelc(k,i)=gelc(k,i)
      &           +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
-     &           *fac_shield(i)*fac_shield(j)   
+     &           *fac_shield(i)**2*fac_shield(j)**2   
             gelc(k,j)=gelc(k,j)
      &           +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
-     &           *fac_shield(i)*fac_shield(j)
+     &           *fac_shield(i)**2*fac_shield(j)**2
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
           enddo
@@ -4151,7 +4166,7 @@ C Contribution to the local-electrostatic energy coming from the i-j pair
      &     +a33*muij(4)
 c          write (iout,*) 'i',i,' j',j,itype(i),itype(j),
 c     &                     ' eel_loc_ij',eel_loc_ij
-c          write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
+C          write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
 C Calculate patrial derivative for theta angle
 #ifdef NEWCORR
          geel_loc_ij=a22*gmuij1(1)
@@ -7289,6 +7304,7 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.SHIELD'
       double precision gx(3),gx1(3)
       logical lprn
       lprn=.false.
@@ -7707,6 +7723,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
       include 'COMMON.CONTACTS'
       include 'COMMON.CHAIN'
       include 'COMMON.CONTROL'
+      include 'COMMON.SHIELD'
       double precision gx(3),gx1(3)
       integer num_cont_hb_old(maxres)
       logical lprn,ldone
@@ -8009,6 +8026,7 @@ cd               write (iout,*) "grij_hb_cont i1",grij_hb_cont(:,jj,i1)
                 call calc_eello(i,jp,i+1,jp1,jj,kk)
                 if (wcorr4.gt.0.0d0) 
      &            ecorr=ecorr+eello4(i,jp,i+1,jp1,jj,kk)
+CC     &            *fac_shield(i)**2*fac_shield(j)**2
                   if (energy_dec.and.wcorr4.gt.0.0d0) 
      1                 write (iout,'(a6,4i5,0pf7.3)')
      2                'ecorr4',i,j,i+1,j1,eello4(i,jp,i+1,jp1,jj,kk)
@@ -8128,6 +8146,7 @@ c------------------------------------------------------------------------------
       include 'COMMON.DERIV'
       include 'COMMON.INTERACT'
       include 'COMMON.CONTACTS'
+      include 'COMMON.SHIELD'
       double precision gx(3),gx1(3)
       logical lprn
       lprn=.false.
@@ -8138,7 +8157,8 @@ c------------------------------------------------------------------------------
       ees0mij=ees0m(jj,i)
       ees0mkl=ees0m(kk,k)
       ekont=eij*ekl
-      ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+      ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)*
+     & fac_shield(i)**2*fac_shield(j)**2
 cd    ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
 C Following 4 lines for diagnostics.
 cd    ees0pkl=0.0D0
@@ -10764,6 +10784,10 @@ C now costhet_grad
      &                    /VSolvSphere_div
 C now the gradient...
 C grad_shield is gradient of Calfa for peptide groups
+C      write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist,
+C     &               costhet,cosphi
+C       write(iout,*) "cosphi_compon",i,k,pep_side0pept_group,
+C     & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k)
       do j=1,3
       grad_shield(j,i)=grad_shield(j,i)
 C gradient po skalowaniu
UNRESPACK 3.x