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
-  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

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
+!        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)
@@ -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)
+!        print *,"before",i,i+3, gshieldc_t4(2,i+3),gshieldc_t4(2,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
-
+      
       j=i+3
+!      if (j.ne.20) return
+!      print *,i,j,gshieldc_t4(2,j),gshieldc_t4(2,j+1)
 !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)
+!           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)
@@ -4635,16 +4640,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)
@@ -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)
+!           print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1)
            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))
+!        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)
@@ -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
+!        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
@@ -10567,7 +10582,7 @@
                      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)
-
+!                      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), &
@@ -10871,7 +10890,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
@@ -10893,6 +10913,7 @@
         write (iout,*) i,gloc(i,icg)
       enddo
 #endif
+!#undef DEBUG
 #ifdef MPI
       if (nfgtasks.gt.1) then
         do j=1,3
@@ -11054,7 +11075,7 @@
         endif
       endif
       endif
-!el#define DEBUG
+!#define DEBUG
 #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
-!el#undef DEBUG
+!#undef DEBUG
 #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
+          call zerograd
           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
+          call zerograd
           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
+        call zerograd
         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
+           call zerograd
             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
+            call zerograd
             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
+            call zerograd
             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
-            call etotal(energia1)
+           call zerograd
+           call etotal(energia1)
             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)
+!              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
@@ -12071,6 +12101,7 @@
         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
@@ -12114,6 +12145,7 @@
 #endif
 !          call int_from_cart1(.false.)
           if (.not.split_ene) then
+           call zerograd
             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 zerograd
             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
+            call zerograd
             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
+            call zerograd
             call etotal(energia1)
             etot2=energia1(0)
           ggg(j+3)=(etot1-etot2)/(2*aincr)
@@ -16212,6 +16247,7 @@
       call sum_gradient
 #ifdef TIMING
 #endif
+!#define DEBUG
 !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
+!#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
@@ -16246,6 +16283,7 @@
 #endif
 !     call checkintcartgrad
 !     write(iout,*) 'calling int_to_cart'
+!#define DEBUG
 #ifdef DEBUG
       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 #endif
@@ -16277,6 +16315,7 @@
                 (gxcart(j,i),j=1,3)
             enddo
 #endif
+!#undef DEBUG
 #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)
-              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
@@ -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)
-               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                                                                                                         
@@ -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)
+!#define DEBUG
 #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
+!#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)
@@ -17048,6 +17094,7 @@
 #endif
             endif
 #endif
+!#define DEBUG
 #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
+!#undef DEBUG
             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
-!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
@@ -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
+!      print *,"SORT",short,long,sinthet,costhet
 !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)) &
-     &                    /VSolvSphere_div
+                         /VSolvSphere_div
 !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
-
+!       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

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       
+!#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)
+!#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)
+!#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
@@ -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)
+!#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                                                                                                                                                   

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
+      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
@@ -3891,22 +3916,22 @@
        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      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
-      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 
-      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