795 |
|
f4 = 0.4_dp*alpha2*f3*r2 |
796 |
|
endif |
797 |
|
ri5damp = f1 + f3*one_third |
798 |
< |
ri7damp = ri5damp + f4 |
798 |
> |
ri7damp = ri5damp + f4*one_third |
799 |
|
|
800 |
|
ri2 = riji * riji |
801 |
|
ri3 = ri2 * riji |
917 |
|
endif |
918 |
|
|
919 |
|
if (j_is_Dipole) then |
920 |
< |
!!$ if (screeningMethod .eq. DAMPED) then |
921 |
< |
!!$ ! assemble the damping variables |
922 |
< |
!!$ call lookupUniformSpline1d(f0spline, rij, f0, df0) |
923 |
< |
!!$ f1 = -rij * df0 + f0 |
924 |
< |
!!$ f2 = -2.0_dp*alpha2*df0 |
925 |
< |
!!$ f3 = f2*r2*rij |
926 |
< |
!!$ f4 = 0.4_dp*alpha2*f3*r2 |
927 |
< |
!!$ endif |
920 |
> |
if (screeningMethod .eq. DAMPED) then |
921 |
> |
! assemble the damping variables |
922 |
> |
call lookupUniformSpline1d(f0spline, rij, f0, df0) |
923 |
> |
f1 = -rij * df0 + f0 |
924 |
> |
f2 = -2.0_dp*alpha2*df0 |
925 |
> |
f3 = f2*r2*rij |
926 |
> |
f4 = 0.4_dp*alpha2*f3*r2 |
927 |
> |
endif |
928 |
|
|
929 |
|
ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3) |
930 |
|
|
981 |
|
scale = 1.0_dp |
982 |
|
endif |
983 |
|
endif |
984 |
– |
|
985 |
– |
sc2 = scale * scale |
984 |
|
|
987 |
– |
pot_term = (ct_ij - 3.0_dp * ct_i * ct_j * sc2) |
988 |
– |
!!$ vterm = pref * ( ri3*pot_term*f1 + (ct_i * ct_j * sc2)*f2 ) |
989 |
– |
vterm = pref * ri3 * pot_term |
990 |
– |
vpair = vpair + vterm |
991 |
– |
epot = epot + sw*vterm |
992 |
– |
|
985 |
|
! precompute variables for convenience (and obfuscation |
986 |
|
! unfortunatly) |
987 |
< |
!!$ ri5damp = f1 + f3*one_third |
988 |
< |
!!$ ri7damp = 5.0_dp*(ri5damp + f4) |
987 |
> |
sc2 = scale * scale |
988 |
> |
ri5damp = f1 + f3*one_third |
989 |
> |
ri7damp = 5.0_dp*(ri5damp + f4*one_third) |
990 |
|
prei3 = sw*pref*ri3 |
991 |
|
prei4 = 3.0_dp*sw*pref*ri4*scale |
992 |
< |
!!$ cti3 = 3.0_dp*ct_i*sc2*ri5damp |
993 |
< |
!!$ ctj3 = 3.0_dp*ct_j*sc2*ri5damp |
994 |
< |
cti3 = 3.0_dp*ct_i*sc2 |
1002 |
< |
ctj3 = 3.0_dp*ct_j*sc2 |
1003 |
< |
ctidotj = ct_i * ct_j * sc2 |
992 |
> |
cti3 = 3.0_dp*ct_i*sc2*ri5damp |
993 |
> |
ctj3 = 3.0_dp*ct_j*sc2*ri5damp |
994 |
> |
ctidotj = ct_i * ct_j * sc2 |
995 |
|
|
996 |
< |
dudx = dudx + prei4 * ( 5.0_dp*ctidotj*xhat - & |
997 |
< |
(ct_i*uz_j(1) + ct_j*uz_i(1) + ct_ij*xhat) ) |
998 |
< |
dudy = dudy + prei4 * ( 5.0_dp*ctidotj*yhat - & |
999 |
< |
(ct_i*uz_j(2) + ct_j*uz_i(2) + ct_ij*yhat) ) |
1000 |
< |
dudz = dudz + prei4 * ( 5.0_dp*ctidotj*zhat - & |
1010 |
< |
(ct_i*uz_j(3) + ct_j*uz_i(3) + ct_ij*zhat) ) |
996 |
> |
! calculate the potential |
997 |
> |
pot_term = (ct_ij*f1 - 3.0_dp*ctidotj*ri5damp) |
998 |
> |
vterm = pref * ri3 * pot_term |
999 |
> |
vpair = vpair + vterm |
1000 |
> |
epot = epot + sw*vterm |
1001 |
|
|
1002 |
< |
duduz_i(1) = duduz_i(1) + prei3 * ( uz_j(1) - ctj3*xhat ) |
1003 |
< |
duduz_i(2) = duduz_i(2) + prei3 * ( uz_j(2) - ctj3*yhat ) |
1004 |
< |
duduz_i(3) = duduz_i(3) + prei3 * ( uz_j(3) - ctj3*zhat ) |
1005 |
< |
|
1006 |
< |
duduz_j(1) = duduz_j(1) + prei3 * ( uz_i(1) - cti3*xhat ) |
1007 |
< |
duduz_j(2) = duduz_j(2) + prei3 * ( uz_i(2) - cti3*yhat ) |
1008 |
< |
duduz_j(3) = duduz_j(3) + prei3 * ( uz_i(3) - cti3*zhat ) |
1002 |
> |
! calculate derivatives for the forces and torques |
1003 |
> |
dudx = dudx + prei4 * ( ctidotj*xhat*ri7damp - & |
1004 |
> |
(ct_i*uz_j(1) + ct_j*uz_i(1) + ct_ij*xhat)*ri5damp ) |
1005 |
> |
dudy = dudy + prei4 * ( ctidotj*yhat*ri7damp - & |
1006 |
> |
(ct_i*uz_j(2) + ct_j*uz_i(2) + ct_ij*yhat)*ri5damp ) |
1007 |
> |
dudz = dudz + prei4 * ( ctidotj*zhat*ri7damp - & |
1008 |
> |
(ct_i*uz_j(3) + ct_j*uz_i(3) + ct_ij*zhat)*ri5damp ) |
1009 |
|
|
1010 |
< |
!!$ dudx = dudx + prei4 * ( ctidotj*xhat*ri7damp - & |
1011 |
< |
!!$ (ct_i*uz_j(1) + ct_j*uz_i(1) + ct_ij*xhat)*ri5damp ) |
1012 |
< |
!!$ dudy = dudy + prei4 * ( ctidotj*yhat*ri7damp - & |
1013 |
< |
!!$ (ct_i*uz_j(2) + ct_j*uz_i(2) + ct_ij*yhat)*ri5damp ) |
1014 |
< |
!!$ dudz = dudz + prei4 * ( ctidotj*zhat*ri7damp - & |
1015 |
< |
!!$ (ct_i*uz_j(3) + ct_j*uz_i(3) + ct_ij*zhat)*ri5damp ) |
1016 |
< |
!!$ |
1027 |
< |
!!$ duduz_i(1) = duduz_i(1) + prei3 * ( uz_j(1)*f1 - ctj3*xhat ) |
1028 |
< |
!!$ duduz_i(2) = duduz_i(2) + prei3 * ( uz_j(2)*f1 - ctj3*yhat ) |
1029 |
< |
!!$ duduz_i(3) = duduz_i(3) + prei3 * ( uz_j(3)*f1 - ctj3*zhat ) |
1030 |
< |
!!$ |
1031 |
< |
!!$ duduz_j(1) = duduz_j(1) + prei3 * ( uz_i(1)*f1 - cti3*xhat ) |
1032 |
< |
!!$ duduz_j(2) = duduz_j(2) + prei3 * ( uz_i(2)*f1 - cti3*yhat ) |
1033 |
< |
!!$ duduz_j(3) = duduz_j(3) + prei3 * ( uz_i(3)*f1 - cti3*zhat ) |
1010 |
> |
duduz_i(1) = duduz_i(1) + prei3 * ( uz_j(1)*f1 - ctj3*xhat ) |
1011 |
> |
duduz_i(2) = duduz_i(2) + prei3 * ( uz_j(2)*f1 - ctj3*yhat ) |
1012 |
> |
duduz_i(3) = duduz_i(3) + prei3 * ( uz_j(3)*f1 - ctj3*zhat ) |
1013 |
> |
|
1014 |
> |
duduz_j(1) = duduz_j(1) + prei3 * ( uz_i(1)*f1 - cti3*xhat ) |
1015 |
> |
duduz_j(2) = duduz_j(2) + prei3 * ( uz_i(2)*f1 - cti3*yhat ) |
1016 |
> |
duduz_j(3) = duduz_j(3) + prei3 * ( uz_i(3)*f1 - cti3*zhat ) |
1017 |
|
|
1018 |
|
endif |
1019 |
|
endif |
1030 |
|
f4 = 0.4_dp*alpha2*f3*r2 |
1031 |
|
endif |
1032 |
|
ri5damp = f1 + f3*one_third |
1033 |
< |
ri7damp = ri5damp + f4 |
1033 |
> |
ri7damp = ri5damp + f4*one_third |
1034 |
|
|
1035 |
|
ri2 = riji * riji |
1036 |
|
ri3 = ri2 * riji |