working gradient for lipid and shielding
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Fri, 7 Jul 2017 08:00:28 +0000 (10:00 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Fri, 7 Jul 2017 08:00:28 +0000 (10:00 +0200)
source/unres/energy.f90

index 014ce41..77c46dd 100644 (file)
 !     if (icall.eq.0) lprn=.false.
 !el      ind=0
       do i=iatsc_s,iatsc_e
-        print *,"I am in EVDW",i
+!C        print *,"I am in EVDW",i
         itypi=iabs(itype(i))
+!        if (i.ne.47) cycle
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
 !el            ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
+!             if (j.ne.78) cycle
 !            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
 !            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,&
     
 
 !d      write(iout,*) 'In EELEC'
-        print *,"IN EELEC"
+!        print *,"IN EELEC"
 !d      do i=1,nloctyp
 !d        write(iout,*) 'Type',i
 !d        write(iout,*) 'B1',B1(:,i)
         time_mat=time_mat+MPI_Wtime()-time01
 #endif
       endif
-       print *, "after set matrices"
+!       print *, "after set matrices"
 !d      do i=1,nres-1
 !d        write (iout,*) 'i=',i
 !d        do k=1,3
 !
 
 
-        print *,"before iturn3 loop"
+!        print *,"before iturn3 loop"
       do i=iturn3_start,iturn3_end
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
         .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
          sslipi=0.0d0
          ssgradlipi=0.0
        endif 
-       print *,i,sslipi,ssgradlipi
+!       print *,i,sslipi,ssgradlipi
        call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
         num_cont_hb(i)=num_conti
 !
 ! Radial derivatives. First process both termini of the fragment (i,j)
 !
-          ggg(1)=facel*xj+sss_ele_grad*rmij*eesij*xj
-          ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj
-          ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj
+          ggg(1)=facel*xj+sss_ele_grad*rmij*eesij*xj* &
+          ((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+          ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj* & 
+           ((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+          ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj* &
+            ((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. &
           (shield_mode.gt.0)) then
             ggg(l)=(agg(l,1)*muij(1)+ &
                 agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))&
             *sss_ele_cut &
-             +eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+          *fac_shield(i)*fac_shield(j) &
+          *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) &
+             +eel_loc_ij*sss_ele_grad*rmij*xtemp(l) 
+
 
             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
           enddo
             gel_loc_long(3,j)=gel_loc_long(3,j)+ &
           ssgradlipj*eel_loc_ij/2.0d0*lipscale/  &
-          ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+          ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
 
             gel_loc_long(3,i)=gel_loc_long(3,i)+ &
           ssgradlipi*eel_loc_ij/2.0d0*lipscale/  &
-          ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+          ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
 
 !grad          do k=i+1,j2
 !grad            do l=1,3
           do l=1,3
             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ &
                 aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))&
-            *sss_ele_cut
+            *sss_ele_cut &
+          *fac_shield(i)*fac_shield(j) &
+          *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 !+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ &
                 aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3) &
             +aggi1(l,4)*muij(4))&
-            *sss_ele_cut
+            *sss_ele_cut &
+          *fac_shield(i)*fac_shield(j) &
+          *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 !+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ &
                 aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))&
-            *sss_ele_cut
+            *sss_ele_cut &
+          *fac_shield(i)*fac_shield(j) &
+          *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 !+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ &
                 aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3) &
             +aggj1(l,4)*muij(4))&
-            *sss_ele_cut
+            *sss_ele_cut &
+          *fac_shield(i)*fac_shield(j) &
+          *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 !+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
           enddo
           ENDIF
         call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1))
         call transpose2(auxmat2(1,1),auxmat3(1,1))
         call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
-        gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
+        gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))&
+          *fac_shield(i)*fac_shield(j)        &
+          *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 ! Derivatives in gamma(i+1)
         call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
         call transpose2(auxmat2(1,1),auxmat3(1,1))
             call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
             s3=0.5d0*(pizda(1,1)+pizda(2,2))
             ggg(l)=-(s1+s2+s3)
-            gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)
+            gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)&
+       *fac_shield(i)*fac_shield(j)  &
+       *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
           enddo
         endif
 ! Remaining derivatives of this turn contribution
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
         time00=MPI_Wtime()
-        call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres,&
+        call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*nres+3,&
           MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres,&
           MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
       do k=1,3
         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss_ele_cut
 !C      print *,'gg',k,gg(k)
-      enddo 
+       enddo 
+!       print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut
 !      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
         gvdwx(k,i)=gvdwx(k,i)-gg(k) +gg_lipi(k)&
@@ -15675,7 +15701,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 #ifdef DEBUG
       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 #endif
-      do i=1,nct
+      do i=0,nct
         do j=1,3
           gcart(j,i)=gradc(j,i,icg)
           gxcart(j,i)=gradx(j,i,icg)
@@ -15703,7 +15729,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 #ifdef DEBUG
       write (iout,*) "CARGRAD"
 #endif
-      do i=nres,1,-1
+      do i=nres,0,-1
         do j=1,3
           gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
 !          gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
@@ -15817,6 +15843,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           gscloc(j,i)=0.0d0
           gsclocx(j,i)=0.0d0
           gliptran(j,i)=0.0d0
+          gliptranx(j,i)=0.0d0
+          gliptranc(j,i)=0.0d0
           gshieldx(j,i)=0.0d0
           gshieldc(j,i)=0.0d0
           gshieldc_loc(j,i)=0.0d0
@@ -18050,7 +18078,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       enddo
       fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
      
-      write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i)
+!C      write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i)
       enddo
       return
       end subroutine set_shield_fac2