ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2438 by chrisfen, Tue Nov 15 19:04:02 2005 UTC vs.
Revision 2439 by chrisfen, Tue Nov 15 19:42:22 2005 UTC

# Line 690 | Line 690 | contains
690         if (j_is_Dipole) then
691  
692            pref = pre12 * q_i * mu_j
693 –
694 – !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
695 – !!$             ri2 = riji * riji
696 – !!$             ri3 = ri2 * riji
697 – !!$
698 – !!$             pref = pre12 * q_i * mu_j
699 – !!$             vterm = - pref * ct_j * (ri2 - rcuti2)
700 – !!$             vpair = vpair + vterm
701 – !!$             epot = epot + sw*vterm
702 – !!$            
703 – !!$             !! this has a + sign in the () because the rij vector is
704 – !!$             !! r_j - r_i and the charge-dipole potential takes the origin
705 – !!$             !! as the point dipole, which is atom j in this case.
706 – !!$            
707 – !!$             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
708 – !!$                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
709 – !!$             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
710 – !!$                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
711 – !!$             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
712 – !!$                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
713 – !!$            
714 – !!$             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
715 – !!$             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
716 – !!$             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
717 – !!$
718 – !!$          elseif (summationMethod .eq. REACTION_FIELD) then
693  
694            if (summationMethod .eq. REACTION_FIELD) then
695               ri2 = riji * riji
696               ri3 = ri2 * riji
697      
724 –             pref = pre12 * q_i * mu_j
698               vterm = - pref * ct_j * ( ri2 - preRF2*rij )
699               vpair = vpair + vterm
700               epot = epot + sw*vterm
# Line 754 | Line 727 | contains
727               ri3 = ri2 * ri
728               sc2 = scale * scale
729  
757 –             pref = pre12 * q_i * mu_j
730               vterm = - pref * ct_j * ri2 * scale
731               vpair = vpair + vterm
732               epot = epot + sw*vterm
# Line 782 | Line 754 | contains
754            cy2 = cy_j * cy_j
755            cz2 = cz_j * cz_j
756  
757 < !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
758 < !!$             pref =  pre14 * q_i / 3.0_dp
759 < !!$             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
760 < !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
761 < !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
762 < !!$             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
791 < !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
792 < !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
793 < !!$             vpair = vpair + ( vterm1 - vterm2 )
794 < !!$             epot = epot + sw*( vterm1 - vterm2 )
795 < !!$            
796 < !!$             dudx = dudx - (5.0_dp * &
797 < !!$                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
798 < !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
799 < !!$                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
800 < !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
801 < !!$                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
802 < !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
803 < !!$                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
804 < !!$             dudy = dudy - (5.0_dp * &
805 < !!$                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
806 < !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
807 < !!$                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
808 < !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
809 < !!$                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
810 < !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
811 < !!$                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
812 < !!$             dudz = dudz - (5.0_dp * &
813 < !!$                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
814 < !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
815 < !!$                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
816 < !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
817 < !!$                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
818 < !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
819 < !!$                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
820 < !!$            
821 < !!$             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
822 < !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
823 < !!$             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
824 < !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
825 < !!$             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
826 < !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
827 < !!$            
828 < !!$             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
829 < !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
830 < !!$             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
831 < !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
832 < !!$             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
833 < !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
834 < !!$            
835 < !!$             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
836 < !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
837 < !!$             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
838 < !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
839 < !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
840 < !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
841 < !!$        
842 < !!$          else
843 <             pref =  pre14 * q_i / 3.0_dp
844 <             vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
845 <                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
846 <                  qzz_j * (3.0_dp*cz2 - 1.0_dp))
847 <             vpair = vpair + vterm
848 <             epot = epot + sw*vterm
849 <            
850 <             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
851 <                  qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
852 <                  qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
853 <                  qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
854 <             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
855 <                  qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
856 <                  qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
857 <                  qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
858 <             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
859 <                  qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
860 <                  qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
861 <                  qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
862 <            
863 <             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
864 <             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
865 <             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
866 <            
867 <             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
868 <             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
869 <             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
870 <            
871 <             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
872 <             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
873 <             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
757 >          pref =  pre14 * q_i / 3.0_dp
758 >          vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
759 >               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
760 >               qzz_j * (3.0_dp*cz2 - 1.0_dp))
761 >          vpair = vpair + vterm
762 >          epot = epot + sw*vterm
763            
764 < !!$          endif
764 >          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
765 >               qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
766 >               qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
767 >               qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
768 >          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
769 >               qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
770 >               qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
771 >               qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
772 >          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
773 >               qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
774 >               qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
775 >               qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
776 >          
777 >          dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
778 >          dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
779 >          dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
780 >          
781 >          duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
782 >          duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
783 >          duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
784 >          
785 >          duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
786 >          duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
787 >          duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
788 >          
789         endif
790      endif
791 <
791 >    
792      if (i_is_Dipole) then
793  
794         if (j_is_Charge) then
795            
796 +          pref = pre12 * q_j * mu_i
797 +          
798            if (summationMethod .eq. SHIFTED_POTENTIAL) then
799               ri2 = riji * riji
800               ri3 = ri2 * riji
801              
887 –             pref = pre12 * q_j * mu_i
802               pot_term = ri2 - rcuti2
803               vterm = pref * ct_i * pot_term
804               vpair = vpair + vterm
# Line 902 | Line 816 | contains
816               ri2 = riji * riji
817               ri3 = ri2 * riji
818  
905 –             pref = pre12 * q_j * mu_i
819               pot_term = ri2 - rcuti2 + 2.0d0*rcuti3*( rij - defaultCutoff )
820               vterm = pref * ct_i * pot_term
821               vpair = vpair + vterm
# Line 920 | Line 833 | contains
833               ri2 = riji * riji
834               ri3 = ri2 * riji
835  
923 –             pref = pre12 * q_j * mu_i
836               vterm = pref * ct_i * ( ri2 - preRF2*rij )
837               vpair = vpair + vterm
838               epot = epot + sw*vterm
# Line 950 | Line 862 | contains
862               ri3 = ri2 * ri
863               sc2 = scale * scale
864  
953 –             pref = pre12 * q_j * mu_i
865               vterm = pref * ct_i * ri2 * scale
866               vpair = vpair + vterm
867               epot = epot + sw*vterm
# Line 974 | Line 885 | contains
885            
886            pref = pre22 * mu_i * mu_j
887  
977 – !!$          if (summationMethod .eq. SHIFTED_POTENTIAL) then
978 – !!$             a0 = ct_ij - 3.0d0 * ct_i * ct_j
979 – !!$             pot_term = ri3 - rcuti3
980 – !!$            
981 – !!$             vterm = pref*pot_term*a0
982 – !!$             vpair = vpair + vterm
983 – !!$             epot = epot + sw*vterm
984 – !!$            
985 – !!$             a1 = 5.0d0 * ct_i * ct_j - ct_ij
986 – !!$            
987 – !!$             dudx = dudx + sw*pref*3.0d0*ri4 &
988 – !!$                  * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
989 – !!$             dudy = dudy + sw*pref*3.0d0*ri4 &
990 – !!$                  * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
991 – !!$             dudz = dudz + sw*pref*3.0d0*ri4 &
992 – !!$                  * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
993 – !!$            
994 – !!$             duduz_i(1) = duduz_i(1) + sw*pref*( pot_term &
995 – !!$                  * (uz_j(1) - 3.0d0*ct_j*xhat) )
996 – !!$             duduz_i(2) = duduz_i(2) + sw*pref*( pot_term &
997 – !!$                  * (uz_j(2) - 3.0d0*ct_j*yhat) )
998 – !!$             duduz_i(3) = duduz_i(3) + sw*pref*( pot_term &
999 – !!$                  * (uz_j(3) - 3.0d0*ct_j*zhat) )
1000 – !!$             duduz_j(1) = duduz_j(1) + sw*pref*( pot_term &
1001 – !!$                  * (uz_i(1) - 3.0d0*ct_i*xhat) )
1002 – !!$             duduz_j(2) = duduz_j(2) + sw*pref*( pot_term &
1003 – !!$                  * (uz_i(2) - 3.0d0*ct_i*yhat) )
1004 – !!$             duduz_j(3) = duduz_j(3) + sw*pref*( pot_term &
1005 – !!$                  * (uz_i(3) - 3.0d0*ct_i*zhat) )
1006 – !!$
1007 – !!$          elseif (summationMethod .eq. SHIFTED_FORCE) then
1008 – !!$             a0 = ct_ij - 3.0d0 * ct_i * ct_j
1009 – !!$             pot_term = ri3 - rcuti3 + 3.0d0*rcuti4*( rij - defaultCutoff )
1010 – !!$            
1011 – !!$             vterm = pref*pot_term*a0
1012 – !!$             vpair = vpair + vterm
1013 – !!$             epot = epot + sw*vterm
1014 – !!$            
1015 – !!$             a1 = 5.0d0 * ct_i * ct_j - ct_ij
1016 – !!$            
1017 – !!$             dudx = dudx + sw*pref*3.0d0*( ri4 - rcuti4 ) &
1018 – !!$                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1019 – !!$             dudy = dudy + sw*pref*3.0d0*( ri4 - rcuti4 ) &
1020 – !!$                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1021 – !!$             dudz = dudz + sw*pref*3.0d0*( ri4 - rcuti4 ) &
1022 – !!$                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1023 – !!$            
1024 – !!$             duduz_i(1) = duduz_i(1) + sw*pref*( pot_term &
1025 – !!$                  * (uz_j(1) - 3.0d0*ct_j*xhat) )
1026 – !!$             duduz_i(2) = duduz_i(2) + sw*pref*( pot_term &
1027 – !!$                  * (uz_j(2) - 3.0d0*ct_j*yhat) )
1028 – !!$             duduz_i(3) = duduz_i(3) + sw*pref*( pot_term &
1029 – !!$                  * (uz_j(3) - 3.0d0*ct_j*zhat) )
1030 – !!$             duduz_j(1) = duduz_j(1) + sw*pref*( pot_term &
1031 – !!$                  * (uz_i(1) - 3.0d0*ct_i*xhat) )
1032 – !!$             duduz_j(2) = duduz_j(2) + sw*pref*( pot_term &
1033 – !!$                  * (uz_i(2) - 3.0d0*ct_i*yhat) )
1034 – !!$             duduz_j(3) = duduz_j(3) + sw*pref*( pot_term &
1035 – !!$                  * (uz_i(3) - 3.0d0*ct_i*zhat) )
1036 – !!$            
1037 – !!$          elseif (summationMethod .eq. REACTION_FIELD) then
888            if (summationMethod .eq. REACTION_FIELD) then
889               vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
890                    preRF2*ct_ij )
# Line 1117 | Line 967 | contains
967  
968      if (i_is_Quadrupole) then
969         if (j_is_Charge) then
1120 –
970            ri2 = riji * riji
971            ri3 = ri2 * riji
972            ri4 = ri2 * ri2
# Line 1125 | Line 974 | contains
974            cy2 = cy_i * cy_i
975            cz2 = cz_i * cz_i
976  
977 < !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
978 < !!$             pref = pre14 * q_j / 3.0_dp
979 < !!$             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
980 < !!$                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
981 < !!$                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
982 < !!$             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
983 < !!$                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
984 < !!$                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
985 < !!$             vpair = vpair + ( vterm1 - vterm2 )
986 < !!$             epot = epot + sw*( vterm1 - vterm2 )
987 < !!$            
988 < !!$             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
989 < !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
990 < !!$                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
991 < !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
992 < !!$                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
993 < !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
994 < !!$                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
995 < !!$             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
996 < !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
997 < !!$                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
998 < !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
999 < !!$                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1000 < !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1001 < !!$                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1002 < !!$             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1003 < !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1004 < !!$                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1005 < !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1006 < !!$                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1007 < !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1008 < !!$                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1160 < !!$            
1161 < !!$             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1162 < !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1163 < !!$             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1164 < !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1165 < !!$             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1166 < !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1167 < !!$            
1168 < !!$             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1169 < !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1170 < !!$             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1171 < !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1172 < !!$             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1173 < !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1174 < !!$            
1175 < !!$             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1176 < !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1177 < !!$             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1178 < !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1179 < !!$             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1180 < !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1181 < !!$
1182 < !!$          else
1183 <             pref = pre14 * q_j / 3.0_dp
1184 <             vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1185 <                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1186 <                  qzz_i * (3.0_dp*cz2 - 1.0_dp))
1187 <             vpair = vpair + vterm
1188 <             epot = epot + sw*vterm
1189 <            
1190 <             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1191 <                  qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1192 <                  qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1193 <                  qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1194 <             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1195 <                  qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1196 <                  qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1197 <                  qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1198 <             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1199 <                  qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1200 <                  qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1201 <                  qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1202 <            
1203 <             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1204 <             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1205 <             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1206 <            
1207 <             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1208 <             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1209 <             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1210 <            
1211 <             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1212 <             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1213 <             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1214 < !!$          endif
977 >          pref = pre14 * q_j / 3.0_dp
978 >          vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
979 >               qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
980 >               qzz_i * (3.0_dp*cz2 - 1.0_dp))
981 >          vpair = vpair + vterm
982 >          epot = epot + sw*vterm
983 >          
984 >          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
985 >               qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
986 >               qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
987 >               qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
988 >          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
989 >               qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
990 >               qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
991 >               qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
992 >          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
993 >               qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
994 >               qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
995 >               qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
996 >          
997 >          dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
998 >          dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
999 >          dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1000 >          
1001 >          duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1002 >          duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1003 >          duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1004 >          
1005 >          duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1006 >          duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1007 >          duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1008 >
1009         endif
1010      endif
1011  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines