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