| 1 |
|
module dipole_dipole |
| 2 |
|
|
| 3 |
< |
use simulation |
| 3 |
> |
use force_globals |
| 4 |
|
use definitions |
| 5 |
|
use atype_module |
| 6 |
|
use vector_class |
| 8 |
|
use mpiSimulation |
| 9 |
|
#endif |
| 10 |
|
implicit none |
| 11 |
< |
|
| 12 |
< |
contains |
| 13 |
< |
|
| 14 |
< |
subroutine do_dipole_pair(atom1, atom2, d, rij, pot, u_l, f, t, & |
| 11 |
> |
|
| 12 |
> |
real(kind=dp), save :: rrf |
| 13 |
> |
real(kind=dp), save :: rt |
| 14 |
> |
real(kind=dp), save :: rrfsq |
| 15 |
> |
real(kind=dp), save :: pre |
| 16 |
> |
logical, save :: dipole_initialized = .false. |
| 17 |
> |
|
| 18 |
> |
contains |
| 19 |
> |
|
| 20 |
> |
subroutine initialize_dipole(this_rrf, this_rt) |
| 21 |
> |
real(kind=dp), intent(in) :: this_rrf, this_rt |
| 22 |
> |
rrf = this_rrf |
| 23 |
> |
rt = this_rt |
| 24 |
> |
rrfsq = rrf * rrf |
| 25 |
> |
! pre converts from mu in units of debye to kcal/mol |
| 26 |
> |
pre = 14.38362d0 |
| 27 |
> |
|
| 28 |
> |
dipole_initialized = .true. |
| 29 |
> |
|
| 30 |
> |
return |
| 31 |
> |
end subroutine initialize_dipole |
| 32 |
> |
|
| 33 |
> |
|
| 34 |
> |
subroutine do_dipole_pair(atom1, atom2, d, rij, r2, pot, u_l, f, t, & |
| 35 |
|
do_pot, do_stress) |
| 36 |
|
|
| 37 |
|
logical :: do_pot, do_stress |
| 38 |
|
|
| 39 |
|
integer atom1, atom2, me1, me2 |
| 40 |
< |
double precision rij, mu1, mu2 |
| 41 |
< |
double precision dfact1, dfact2, dip2, r2, r3, r5, pre |
| 42 |
< |
double precision dudx, dudy, dudz, dudu1x, dudu1y, dudu1z |
| 43 |
< |
double precision dudu2x, dudu2y, dudu2z, rdotu1, rdotu2, u1dotu2 |
| 44 |
< |
double precision taper, dtdr, vterm |
| 40 |
> |
real(kind=dp) :: rij, mu1, mu2 |
| 41 |
> |
real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5, pre |
| 42 |
> |
real(kind=dp) :: dudx, dudy, dudz, dudu1x, dudu1y, dudu1z |
| 43 |
> |
real(kind=dp) :: dudu2x, dudu2y, dudu2z, rdotu1, rdotu2, u1dotu2 |
| 44 |
> |
real(kind=dp) :: taper, dtdr, vterm |
| 45 |
|
|
| 46 |
< |
real( kind = dp ) :: pot, rt, rrf |
| 46 |
> |
real( kind = dp ) :: pot |
| 47 |
|
real( kind = dp ), dimension(3) :: d |
| 48 |
|
real( kind = dp ), dimension(3,getNlocal()) :: u_l |
| 49 |
|
real( kind = dp ), dimension(3,getNlocal()) :: f |
| 52 |
|
real (kind = dp), dimension(3) :: ul1 |
| 53 |
|
real (kind = dp), dimension(3) :: ul2 |
| 54 |
|
|
| 55 |
+ |
if (.not.dipole_initialized) then |
| 56 |
+ |
write(default_error,*) 'Dipole-dipole not initialized!' |
| 57 |
+ |
return |
| 58 |
+ |
endif |
| 59 |
|
|
| 60 |
|
#ifdef IS_MPI |
| 61 |
|
me1 = atid_Row(atom1) |
| 81 |
|
|
| 82 |
|
call getElementProperty(atypes, me1, "dipole_moment", mu1) |
| 83 |
|
call getElementProperty(atypes, me2, "dipole_moment", mu2) |
| 60 |
– |
|
| 61 |
– |
rrf = getRrf() |
| 62 |
– |
rt = getRt() |
| 63 |
– |
|
| 64 |
– |
! pre converts from mu in units of debye to kcal/mol |
| 65 |
– |
pre = 14.38362d0 |
| 84 |
|
|
| 85 |
|
if (rij.le.rrf) then |
| 86 |
|
|
| 92 |
|
dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3) |
| 93 |
|
endif |
| 94 |
|
|
| 77 |
– |
r2 = rij*rij |
| 95 |
|
r3 = r2*rij |
| 96 |
|
r5 = r3*r2 |
| 97 |
|
|
