changes in shielding

diff --git a/PARAM/scparm.adam_30_10_2013 b/PARAM/scparm.adam_30_10_2013
index 82dacd1..d5e645c 100644 (file)
--- a/PARAM/scparm.adam_30_10_2013
+++ b/PARAM/scparm.adam_30_10_2013
@@ -350,7 +350,7 @@
   0.0347  4.9555 -0.1115  0.4550 -0.8819  0.7342  0.0000  0.0000  0.0000  0.0000  3.9312  7.0874 -0.1922  0.1141  1  0.0000  0.0000  0.0000  0.0000  1.6344  0.0000  0.0000  0.0000  0.5564  0.0531  0.0000  5.0000  0.2736  1.5702  0.0000  0.0000  0.0000  1.1690  0.0000  7.6197  0.1701  1.0347 13.6260  0.0000  1.0000  1.0000
   0.0347  4.9555 -0.1115  0.4550 -0.8819  0.7342  0.0000  0.0000  0.0000  0.0000  3.9312  7.0874 -0.1922  0.1141  1  0.0000  0.0000  0.0000  0.0000  1.6344  0.0000  0.0000  0.0000  0.5564  0.0531  0.0000  5.0000  0.2736  1.5702  0.0000  0.0000  0.0000  1.1690  0.0000  7.6197  0.1701  1.0347 13.6260  0.0000  1.0000  1.0000
   0.4409  7.0124 -0.7792 -0.9026 -0.9824  0.9981  0.0000  0.0000  0.0000  0.0000  4.2163  6.7741 -0.9979 -0.6323  1  0.0000  0.0000  0.0000  0.0000  1.2736  0.0000  0.0000  0.0000  6.9619 -0.0105  0.0000  5.0000  0.3624  1.4965  0.0000  0.0000  0.0000  1.5081  0.0000  6.4201  0.1437  0.1181 13.7010  0.0000  1.0000  1.0000
   0.0347  4.9555 -0.1115  0.4550 -0.8819  0.7342  0.0000  0.0000  0.0000  0.0000  3.9312  7.0874 -0.1922  0.1141  1  0.0000  0.0000  0.0000  0.0000  1.6344  0.0000  0.0000  0.0000  0.5564  0.0531  0.0000  5.0000  0.2736  1.5702  0.0000  0.0000  0.0000  1.1690  0.0000  7.6197  0.1701  1.0347 13.6260  0.0000  1.0000  1.0000
   0.0347  4.9555 -0.1115  0.4550 -0.8819  0.7342  0.0000  0.0000  0.0000  0.0000  3.9312  7.0874 -0.1922  0.1141  1  0.0000  0.0000  0.0000  0.0000  1.6344  0.0000  0.0000  0.0000  0.5564  0.0531  0.0000  5.0000  0.2736  1.5702  0.0000  0.0000  0.0000  1.1690  0.0000  7.6197  0.1701  1.0347 13.6260  0.0000  1.0000  1.0000
   0.4409  7.0124 -0.7792 -0.9026 -0.9824  0.9981  0.0000  0.0000  0.0000  0.0000  4.2163  6.7741 -0.9979 -0.6323  1  0.0000  0.0000  0.0000  0.0000  1.2736  0.0000  0.0000  0.0000  6.9619 -0.0105  0.0000  5.0000  0.3624  1.4965  0.0000  0.0000  0.0000  1.5081  0.0000  6.4201  0.1437  0.1181 13.7010  0.0000  1.0000  1.0000

-  1.2808  3.6212  0.2182  0.8105  0.9442 -0.8842  0.0000  0.0000  0.0000  0.0000  1.8965  3.4196  0.1846  0.7056  1  0.0000  0.0000  0.0000  0.0000  1.2281  0.0000  0.0000  0.0000 -0.0631 -0.2997  0.0000  5.0000  1.1189  1.3499  0.0000  0.0000  0.0000  1.3003  0.0000 13.2550 11.0080  0.0349 13.4380  0.0000  1.0000  1.0000  !OThr 1/8 epsilon i +1A sigma
+  1.2808  2.6212  0.2182  0.8105  0.9442 -0.8842  0.0000  0.0000  0.0000  0.0000  1.8965  3.4196  0.1846  0.7056  1  0.0000  0.0000  0.0000  0.0000  1.2281  0.0000  0.0000  0.0000 -0.0631 -0.2997  0.0000  5.0000  1.1189  1.3499  0.0000  0.0000  0.0000  1.3003  0.0000 13.2550 11.0080  0.0349 13.4380  0.0000  1.0000  1.0000  !OThr 1/8 epsilon i 

   0.0354  4.5650  0.4668  0.4668  0.8709  0.8709  0.0000  0.0000  0.0000  0.0000  6.2755  6.2755  0.0232  0.0232  1  0.0000  0.0000  0.0000  0.0000  2.0870  0.0000  0.0000  0.0000  0.0538  0.0538  0.0000  5.0000  1.7567  1.7567  0.0000  0.0000  0.0000  0.0000  0.0000  7.6397  0.1745  1.1729 13.6680  0.0000  1.0000  1.0000
   0.1567  4.6094  0.4496  0.1023 -0.9998 -0.1723  0.0000  0.0000  0.0000  0.0000  2.4369  3.9490 -0.8774  0.3148  1  0.0000  0.0000  0.0000  0.0000  2.8745  0.0000  0.0000  0.0000  3.0813  0.2333  0.0000  5.0000  1.8847  2.0403  0.0000  0.0000  0.0000  1.7847  0.0000 12.0140  7.4670  6.1026 13.6260  0.0000  1.0000  1.0000
   0.6738  4.5858 -0.0217  0.3581 -0.9436  0.8502  0.0000  0.0000  0.0000  0.0000  3.1814  5.1114 -0.9830 -0.9396  1  0.0000  0.0000  0.0000  0.0000  3.0907  0.0000  0.0000  0.0000  1.2890  0.1120  0.0000  5.0000  0.4529  1.6272  0.0000  0.0000  0.0000  1.6273  0.0000  7.4321  0.0349  1.5069 13.7270  0.0000  1.0000  1.0000
   0.0354  4.5650  0.4668  0.4668  0.8709  0.8709  0.0000  0.0000  0.0000  0.0000  6.2755  6.2755  0.0232  0.0232  1  0.0000  0.0000  0.0000  0.0000  2.0870  0.0000  0.0000  0.0000  0.0538  0.0538  0.0000  5.0000  1.7567  1.7567  0.0000  0.0000  0.0000  0.0000  0.0000  7.6397  0.1745  1.1729 13.6680  0.0000  1.0000  1.0000
   0.1567  4.6094  0.4496  0.1023 -0.9998 -0.1723  0.0000  0.0000  0.0000  0.0000  2.4369  3.9490 -0.8774  0.3148  1  0.0000  0.0000  0.0000  0.0000  2.8745  0.0000  0.0000  0.0000  3.0813  0.2333  0.0000  5.0000  1.8847  2.0403  0.0000  0.0000  0.0000  1.7847  0.0000 12.0140  7.4670  6.1026 13.6260  0.0000  1.0000  1.0000
   0.6738  4.5858 -0.0217  0.3581 -0.9436  0.8502  0.0000  0.0000  0.0000  0.0000  3.1814  5.1114 -0.9830 -0.9396  1  0.0000  0.0000  0.0000  0.0000  3.0907  0.0000  0.0000  0.0000  1.2890  0.1120  0.0000  5.0000  0.4529  1.6272  0.0000  0.0000  0.0000  1.6273  0.0000  7.4321  0.0349  1.5069 13.7270  0.0000  1.0000  1.0000

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index fa701a0..e4d0f76 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -3159,6 +3159,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
         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)
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)

@@ -3201,6 +3202,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)
         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
 !
         num_cont_hb(i)=num_conti
       enddo   ! i
 !

@@ -4545,8 +4547,10 @@
       integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield
       real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,&
          rlocshield
       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
       j=i+3

+!      if (j.ne.20) return
+!      print *,i,j,gshieldc_t4(2,j),gshieldc_t4(2,j+1)

 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 !
 !               Fourth-order contributions
 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 !
 !               Fourth-order contributions

@@ -4625,6 +4629,7 @@
            iresshield=shield_list(ilist,i)
            do k=1,3
            rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i)
            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)
            gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ &
                    rlocshield &
             +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i)

@@ -4635,16 +4640,17 @@
           do ilist=1,ishield_list(j)
            iresshield=shield_list(ilist,j)
            do k=1,3
           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
            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
 
            enddo
           enddo

-

           do k=1,3
             gshieldc_t4(k,i)=gshieldc_t4(k,i)+  &
                    grad_shield(k,i)*eello_t4/fac_shield(i)
           do k=1,3
             gshieldc_t4(k,i)=gshieldc_t4(k,i)+  &
                    grad_shield(k,i)*eello_t4/fac_shield(i)

@@ -4654,6 +4660,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)
                    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
 
            enddo
            endif
 

@@ -4772,10 +4779,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))
           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)
           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)
 
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)

@@ -4791,10 +4799,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
           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)
           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
         enddo
          gshieldc_t4(3,i)=gshieldc_t4(3,i)+ &
           ssgradlipi*eello_t4/4.0d0*lipscale

@@ -10567,7 +10582,7 @@
                      wscpho*gvdwc_scpho(j,i)+   &
                      wpeppho*gvdwc_peppho(j,i)
 
                      wscpho*gvdwc_scpho(j,i)+   &
                      wpeppho*gvdwc_peppho(j,i)
 

-
+       

 
 
 
 
 
 

@@ -10765,7 +10780,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)
                      +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), &
 !                 if ((i.le.2).and.(i.ge.1))
 !                       print *,gradc(j,i,icg),&
 !                      gradbufc(j,i),welec*gelc(j,i), &

@@ -10871,7 +10890,8 @@
 !              if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i)
 
         enddo
 !              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
 #ifdef DEBUG
       write (iout,*) "gloc before adding corr"
       do i=1,4*nres

@@ -10893,6 +10913,7 @@
         write (iout,*) i,gloc(i,icg)
       enddo
 #endif
         write (iout,*) i,gloc(i,icg)
       enddo
 #endif

+!#undef DEBUG

 #ifdef MPI
       if (nfgtasks.gt.1) then
         do j=1,3
 #ifdef MPI
       if (nfgtasks.gt.1) then
         do j=1,3

@@ -11054,7 +11075,7 @@
         endif
       endif
       endif
         endif
       endif
       endif

-!el#define DEBUG
+!#define DEBUG

 #ifdef DEBUG
       write (iout,*) "gradc gradx gloc"
       do i=1,nres
 #ifdef DEBUG
       write (iout,*) "gradc gradx gloc"
       do i=1,nres

@@ -11062,7 +11083,7 @@
          i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg)
       enddo 
 #endif
          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
 #ifdef TIMING
       time_sumgradient=time_sumgradient+MPI_Wtime()-time01
 #endif

@@ -11744,6 +11765,7 @@
           c(j,k)=c(j,k)+aincr
           c(j,k+nres)=c(j,k+nres)+aincr
           enddo
           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
           call etotal(energia1)
           etot1=energia1(0)
         ggg(j)=(etot1-etot)/aincr

@@ -11756,6 +11778,7 @@
       do j=1,3
         c(j,i+nres)=c(j,i+nres)+aincr
         dc(j,i+nres)=dc(j,i+nres)+aincr
       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
           call etotal(energia1)
           etot1=energia1(0)
         ggg(j+3)=(etot1-etot)/aincr

@@ -11816,6 +11839,7 @@
       write (iout,*) "split_ene ",split_ene
       call flush(iout)
       if (.not.split_ene) then
       write (iout,*) "split_ene ",split_ene
       call flush(iout)
       if (.not.split_ene) then

+        call zerograd

         call etotal(energia)
         etot=energia(0)
         call cartgrad
         call etotal(energia)
         etot=energia(0)
         call cartgrad

@@ -11892,6 +11916,7 @@
           dc(j,i+nres)=c(j,i+nres)-c(j,i)
           call int_from_cart1(.false.)
           if (.not.split_ene) then
           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
             call etotal(energia1)
             etot1=energia1(0)
             write (iout,*) "ij",i,j," etot1",etot1

@@ -11912,6 +11937,7 @@
           dc(j,i+nres)=c(j,i+nres)-c(j,i)
           call int_from_cart1(.false.)
           if (.not.split_ene) then
           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
             call etotal(energia1)
             etot2=energia1(0)
             write (iout,*) "ij",i,j," etot2",etot2

@@ -11943,6 +11969,7 @@
           dc(j,i+nres)=c(j,i+nres)-c(j,i)
           call int_from_cart1(.false.)
           if (.not.split_ene) then
           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
             call etotal(energia1)
             etot1=energia1(0)
           else

@@ -11957,7 +11984,8 @@
           dc(j,i+nres)=c(j,i+nres)-c(j,i)
           call int_from_cart1(.false.)
           if (.not.split_ene) then
           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
             etot2=energia1(0)
           ggg(j+3)=(etot1-etot2)/(2*aincr)
           else

@@ -12050,6 +12078,8 @@
         do i=1,nres
           do j=1,3
             grad_s(j,i)=gcart(j,i)
         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
 !            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

@@ -12071,6 +12101,7 @@
         do i=1,nres
           do j=1,3
             grad_s(j,i)=gcart(j,i)
         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
             grad_s(j+3,i)=gxcart(j,i)
           enddo
         enddo

@@ -12114,6 +12145,7 @@
 #endif
 !          call int_from_cart1(.false.)
           if (.not.split_ene) then
 #endif
 !          call int_from_cart1(.false.)
           if (.not.split_ene) then

+           call zerograd

             call etotal(energia1)
             etot1=energia1(0)
 !            call enerprint(energia1)
             call etotal(energia1)
             etot1=energia1(0)
 !            call enerprint(energia1)

@@ -12131,6 +12163,7 @@
           call chainbuild_cart
 !          call int_from_cart1(.false.)
           if (.not.split_ene) then
           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)
             call etotal(energia1)
             etot2=energia1(0)
           ggg(j)=(etot1-etot2)/(2*aincr)

@@ -12161,6 +12194,7 @@
 !     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
 !          write (iout,*)
           if (.not.split_ene) then
 !     &      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
             call etotal(energia1)
             etot1=energia1(0)
           else

@@ -12183,6 +12217,7 @@
 !          write (iout,*) "dxnormnormsafe",dsqrt(
 !     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
           if (.not.split_ene) then
 !          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)
             call etotal(energia1)
             etot2=energia1(0)
           ggg(j+3)=(etot1-etot2)/(2*aincr)

@@ -16212,6 +16247,7 @@
       call sum_gradient
 #ifdef TIMING
 #endif
       call sum_gradient
 #ifdef TIMING
 #endif

+!#define DEBUG

 !el      write (iout,*) "After sum_gradient"
 #ifdef DEBUG
 !el      write (iout,*) "After sum_gradient"
 !el      write (iout,*) "After sum_gradient"
 #ifdef DEBUG
 !el      write (iout,*) "After sum_gradient"

@@ -16220,6 +16256,7 @@
         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
       enddo
 #endif
         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
 ! If performing constraint dynamics, add the gradients of the constraint energy
       if(usampl.and.totT.gt.eq_time) then
          do i=1,nct

@@ -16246,6 +16283,7 @@
 #endif
 !     call checkintcartgrad
 !     write(iout,*) 'calling int_to_cart'
 #endif
 !     call checkintcartgrad
 !     write(iout,*) 'calling int_to_cart'

+!#define DEBUG

 #ifdef DEBUG
       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 #endif
 #ifdef DEBUG
       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 #endif

@@ -16277,6 +16315,7 @@
                 (gxcart(j,i),j=1,3)
             enddo
 #endif
                 (gxcart(j,i),j=1,3)
             enddo
 #endif

+!#undef DEBUG

 #ifdef CARGRAD
 #ifdef DEBUG
             write (iout,*) "CARGRAD"
 #ifdef CARGRAD
 #ifdef DEBUG
             write (iout,*) "CARGRAD"

@@ -16601,22 +16640,22 @@
               dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+ &
               cost1*dc_norm(j,i-2))/ &
               vbld(i-1)
               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)
               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)
       !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)
               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
             enddo
              endif
             enddo

@@ -16677,15 +16716,20 @@
                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)
                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)
                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)
                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                                                                                                         
                endif
              enddo
             endif                                                                                                         

@@ -17023,6 +17067,7 @@
             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)
             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)
 #ifdef DEBUG
       !d      write (iout,*) "Gather dphi"
       !d      call flush(iout)

@@ -17031,6 +17076,7 @@
             write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
             enddo
 #endif
             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)
             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)

@@ -17048,6 +17094,7 @@
 #endif
             endif
 #endif
 #endif
             endif
 #endif

+!#define DEBUG

 #ifdef DEBUG
             write (iout,*) "dtheta after gather"
             do i=1,nres
 #ifdef DEBUG
             write (iout,*) "dtheta after gather"
             do i=1,nres

@@ -17066,6 +17113,7 @@
             write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3)
             enddo
 #endif
             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
       !-----------------------------------------------------------------------------
             return
             end subroutine intcartderiv
       !-----------------------------------------------------------------------------

@@ -19367,7 +19415,7 @@
       enddo
 !C now sscale fraction
        sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
       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
 !C now sscale
         if (sh_frac_dist.le.0.0) cycle
 !C        print *,ishield_list(i),i

@@ -19401,6 +19449,7 @@
       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
       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
 !C now costhet_grad
 !C       costhet=0.6d0
 !C       sinthet=0.8

@@ -19450,7 +19499,7 @@
        enddo
 !C      print *,sinphi,sinthet
       VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) &
        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
 !C     &                    *wshield
 !C now the gradient...
       do j=1,3

@@ -19472,19 +19521,22 @@
             sinphi/sinthet*costhet*costhet_grad(j)&
            +sinthet/sinphi*cosphi*cosphi_grad_long(j))) &
             )*wshield
             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
        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)
      
       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
       enddo
       return
       end subroutine set_shield_fac2

diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90
index 2287f0a..f41bd30 100644 (file)
--- a/source/unres/geometry.f90
+++ b/source/unres/geometry.f90
@@ -3331,13 +3331,38 @@
          enddo
       endif 
 !  Settind dE/ddnres-1       
          enddo
       endif 
 !  Settind dE/ddnres-1       

+!#define DEBUG
+#ifdef DEBUG
+          j=1
+              write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
+        gloc(2*nres-5,icg),dtheta(j,2,nres)
+
+#endif
+!#undef DEBUG
+

        do j=1,3
         gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
        gloc(2*nres-5,icg)*dtheta(j,2,nres)
        do j=1,3
         gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
        gloc(2*nres-5,icg)*dtheta(j,2,nres)

+!#define DEBUG
+#ifdef DEBUG
+              write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
+        gloc(2*nres-5,icg),dtheta(j,2,nres)
+
+#endif
+!#undef DEBUG

         if(itype(nres-1,1).ne.10) then
           gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
          dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
           domega(j,2,nres-1)
         if(itype(nres-1,1).ne.10) then
           gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
          dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
           domega(j,2,nres-1)

+!#define DEBUG
+#ifdef DEBUG
+              write(iout,*)"in int to cart2",i,gcart(j,nres-1),gloc(ialph(nres-1,1),icg)* &
+          dalpha(j,2,nres-1),gloc(ialph(nres-1,1)+nside,icg), &
+          domega(j,2,nres-1)
+
+#endif
+!#undef DEBUG
+

         endif
         enddo
 !   The side-chain vector derivatives
         endif
         enddo
 !   The side-chain vector derivatives

@@ -3347,6 +3372,12 @@
             do j=1,3   
               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
               +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
             do j=1,3   
               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
               +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)

+!#define DEBUG
+#ifdef DEBUG
+              write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), &
+              gloc(ialph(i,1)+nside,icg),domega(j,3,i)
+#endif
+!#undef DEBUG

             enddo
          endif     
        enddo                                                                                                                                                   
             enddo
          endif     
        enddo                                                                                                                                                   

diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90
index 7534c41..b3492fe 100644 (file)
--- a/source/unres/io_config.f90
+++ b/source/unres/io_config.f90
@@ -2771,6 +2771,31 @@
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
         ' v3ss:',v3ss
       endif
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
         ' v3ss:',v3ss
       endif

+      if (shield_mode.gt.0) then
+      pi=4.0D0*datan(1.0D0)
+!C VSolvSphere the volume of solving sphere
+      print *,pi,"pi"
+!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
+!C there will be no distinction between proline peptide group and normal peptide
+!C group in case of shielding parameters
+      VSolvSphere=4.0/3.0*pi*(4.50d0)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
+      write (iout,*) VSolvSphere,VSolvSphere_div
+!C long axis of side chain 
+      do i=1,ntyp
+      long_r_sidechain(i)=vbldsc0(1,i)
+!      if (scelemode.eq.0) then
+      short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
+      if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
+!      else
+!      short_r_sidechain(i)=sigma(i,i)
+!      endif
+      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
+         sigma0(i) 
+      enddo
+      buff_shield=1.0d0
+      endif
+

       return
   111 write (iout,*) "Error reading bending energy parameters."
       goto 999
       return
   111 write (iout,*) "Error reading bending energy parameters."
       goto 999

@@ -3891,22 +3916,22 @@
        write (iout,'(2a)') diagmeth(kdiag),&
         ' routine used to diagonalize matrices.'
       if (shield_mode.gt.0) then
        write (iout,'(2a)') diagmeth(kdiag),&
         ' routine used to diagonalize matrices.'
       if (shield_mode.gt.0) then

-      pi=3.141592d0
+       pi=4.0D0*datan(1.0D0)

 !C VSolvSphere the volume of solving sphere
 !C VSolvSphere the volume of solving sphere

-!C      print *,pi,"pi"
+      print *,pi,"pi"

 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
 !C there will be no distinction between proline peptide group and normal peptide
 !C group in case of shielding parameters
 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
 !C there will be no distinction between proline peptide group and normal peptide
 !C group in case of shielding parameters

-      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
-      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+      VSolvSphere=4.0/3.0*pi*(4.50d0)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3

       write (iout,*) VSolvSphere,VSolvSphere_div
 !C long axis of side chain 
       write (iout,*) VSolvSphere,VSolvSphere_div
 !C long axis of side chain 

-      do i=1,ntyp
-      long_r_sidechain(i)=vbldsc0(1,i)
-      short_r_sidechain(i)=sigma0(i)
-      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
-         sigma0(i) 
-      enddo
+!      do i=1,ntyp
+!      long_r_sidechain(i)=vbldsc0(1,i)
+!      short_r_sidechain(i)=sigma0(i)
+!      write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
+!         sigma0(i) 
+!      enddo

       buff_shield=1.0d0
       endif
       return
       buff_shield=1.0d0
       endif
       return