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