working gradient for ecorr (wcorrh)
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Thu, 27 Apr 2017 12:13:13 +0000 (14:13 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Thu, 27 Apr 2017 12:13:13 +0000 (14:13 +0200)
source/unres/energy.f90

index d9e28e6..a7e3576 100644 (file)
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
           enddo
- 128      continue
+
           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
               .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 &
               .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
 !          if (energy_dec) write (iout,*) "muij",muij
 !              write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
 
-          eel_loc=eel_loc+eel_loc_ij
+          eel_loc=eel_loc+eel_loc_ij*sss_ele_cut
 ! Partial derivatives in virtual-bond dihedral angles gamma
           if (i.gt.1) &
           gel_loc_loc(i-1)=gel_loc_loc(i-1)+ &
-                  a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) &
-                 +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
+                  (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) &
+                 +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) &
+                 *sss_ele_cut
           gel_loc_loc(j-1)=gel_loc_loc(j-1)+ &
-                  a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) &
-                 +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
+                  (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) &
+                 +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) &
+                 *sss_ele_cut
 ! Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
-          do l=1,3
-            ggg(l)=agg(l,1)*muij(1)+ &
-                agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+!          do l=1,3
+            ggg(1)=(agg(1,1)*muij(1)+ &
+                agg(1,2)*muij(2)+agg(1,3)*muij(3)+agg(1,4)*muij(4)) &
+            *sss_ele_cut+eel_loc_ij*sss_ele_grad*rmij*xj
+            ggg(2)=(agg(2,1)*muij(1)+ &
+                agg(2,2)*muij(2)+agg(2,3)*muij(3)+agg(2,4)*muij(4)) &
+            *sss_ele_cut+eel_loc_ij*sss_ele_grad*rmij*yj
+            ggg(3)=(agg(3,1)*muij(1)+ &
+                agg(3,2)*muij(2)+agg(3,3)*muij(3)+agg(3,4)*muij(4)) &
+            *sss_ele_cut+eel_loc_ij*sss_ele_grad*rmij*zj
+
+           do l=1,3
             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
 !grad            ghalf=0.5d0*ggg(l)
 !grad          enddo
 ! Remaining derivatives of eello
           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)
-            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)
-            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)
-            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)
+            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
+            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
+            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
+            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
           enddo
           ENDIF
 ! Change 12/26/95 to calculate four-body contributions to H-bonding energy
                   ees0mij=0
                 endif
 !               ees0mij=0.0D0
-                ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
-                ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+                ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
+                     *sss_ele_cut
+
+                ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
+                     *sss_ele_cut
+
 ! Diagnostics. Comment out or remove after debugging!
 !               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
 !               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
                 enddo
-                gggp(1)=gggp(1)+ees0pijp*xj
-                gggp(2)=gggp(2)+ees0pijp*yj
-                gggp(3)=gggp(3)+ees0pijp*zj
-                gggm(1)=gggm(1)+ees0mijp*xj
-                gggm(2)=gggm(2)+ees0mijp*yj
-                gggm(3)=gggm(3)+ees0mijp*zj
+                gggp(1)=gggp(1)+ees0pijp*xj &
+                  +ees0p(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+                gggp(2)=gggp(2)+ees0pijp*yj &
+               +ees0p(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+                gggp(3)=gggp(3)+ees0pijp*zj &
+               +ees0p(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
+                gggm(1)=gggm(1)+ees0mijp*xj &
+               +ees0m(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+
+                gggm(2)=gggm(2)+ees0mijp*yj &
+               +ees0m(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+
+                gggm(3)=gggm(3)+ees0mijp*zj &
+               +ees0m(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
 ! Derivatives due to the contact function
                 gacont_hbr(1,num_conti,i)=fprimcont*xj
                 gacont_hbr(2,num_conti,i)=fprimcont*yj
 !grad                  ghalfm=0.5D0*gggm(k)
                   gacontp_hb1(k,num_conti,i)= & !ghalfp+
                     (ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
-                   + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+                   + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+                     *sss_ele_cut
+
                   gacontp_hb2(k,num_conti,i)= & !ghalfp+
                     (ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
-                   + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
-                  gacontp_hb3(k,num_conti,i)=gggp(k)
+                   + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+                     *sss_ele_cut
+
+                  gacontp_hb3(k,num_conti,i)=gggp(k) &
+                     *sss_ele_cut
+
                   gacontm_hb1(k,num_conti,i)= & !ghalfm+
                     (ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
-                   + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+                   + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+                     *sss_ele_cut
+
                   gacontm_hb2(k,num_conti,i)= & !ghalfm+
                     (ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
-                   + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
-                  gacontm_hb3(k,num_conti,i)=gggm(k)
+                   + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) &
+                     *sss_ele_cut
+
+                  gacontm_hb3(k,num_conti,i)=gggm(k) &
+                     *sss_ele_cut
+
                 enddo
 ! Diagnostics. Comment out or remove after debugging!
 !diag           do k=1,3
               enddo
             endif
           endif
+ 128  continue
 !          t_eelecij=t_eelecij+MPI_Wtime()-time00
       return
       end subroutine eelecij