-
Notifications
You must be signed in to change notification settings - Fork 103
/
test_pot_dq.f90
148 lines (121 loc) · 6.3 KB
/
test_pot_dq.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
! test_pot_dq.f90
! dipole-quadrupole potential
MODULE test_pot_module
!------------------------------------------------------------------------------------------------!
! This software was written in 2016/17 !
! by Michael P. Allen <[email protected]>/<[email protected]> !
! and Dominic J. Tildesley <[email protected]> ("the authors"), !
! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), !
! published by Oxford University Press ("the publishers"). !
! !
! LICENCE !
! Creative Commons CC0 Public Domain Dedication. !
! To the extent possible under law, the authors have dedicated all copyright and related !
! and neighboring rights to this software to the PUBLIC domain worldwide. !
! This software is distributed without any warranty. !
! You should have received a copy of the CC0 Public Domain Dedication along with this software. !
! If not, see <http://creativecommons.org/publicdomain/zero/1.0/>. !
! !
! DISCLAIMER !
! The authors and publishers make no warranties about the software, and disclaim liability !
! for all uses of the software, to the fullest extent permitted by applicable law. !
! The authors and publishers do not recommend use of this software for any purpose. !
! It is made freely available, solely to clarify points made in the text. When using or citing !
! the software, you should not imply endorsement by the authors or publishers. !
!------------------------------------------------------------------------------------------------!
USE, INTRINSIC :: iso_fortran_env, ONLY : error_unit
IMPLICIT NONE
PRIVATE
! Public routine
PUBLIC :: force
! Public data
INTEGER, PARAMETER, PUBLIC :: n = 2 ! pair potential
CONTAINS
SUBROUTINE force ( r, e, pot, f, t )
USE maths_module, ONLY : cross_product
IMPLICIT NONE
REAL, DIMENSION(:,:), INTENT(in) :: r, e
REAL, INTENT(out) :: pot
REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: f, t
REAL, DIMENSION(3) :: rij, sij, fij, gi, gj
REAL :: rij_mag, ci, cj, cij
REAL :: vij_dq, vij_qd, dvdrij, dvdci, dvdcj, dvdcij ! Potential derivatives
REAL, PARAMETER :: tol = 1.e-6
INTEGER, PARAMETER :: i = 1, j = 2 ! notation to match appendix
! Routine to demonstrate the calculation of forces and torques from the
! dipole-quadrupole potential
! Written for ease of comparison with the text, rather than efficiency!
! Check dimensions to be sure
IF ( ANY ( SHAPE(r) /= [3,n] ) ) THEN
WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'r shape error', SHAPE(r), 3, n
STOP 'Error in test_pot_dq'
END IF
IF ( ANY ( SHAPE(e) /= [3,n] ) ) THEN
WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'e shape error', SHAPE(e), 3, n
STOP 'Error in test_pot_dq'
END IF
! Check unit vectors
IF ( ABS(SUM(e(:,i)**2)-1.0) > tol .OR. ABS(SUM(e(:,j)**2)-1.0) > tol ) THEN
WRITE ( unit=error_unit, fmt='(a)' ) 'Warning, non-unit vectors'
WRITE ( unit=error_unit, fmt='(4f10.5)' ) e(:,i), SUM(e(:,i)**2)
WRITE ( unit=error_unit, fmt='(4f10.5)' ) e(:,j), SUM(e(:,j)**2)
END IF
rij = r(:,i) - r(:,j)
rij_mag = SQRT ( SUM(rij**2) ) ! Magnitude of separation vector
sij = rij / rij_mag ! Unit vector
ci = DOT_PRODUCT ( e(:,i), sij )
cj = DOT_PRODUCT ( e(:,j), sij )
cij = DOT_PRODUCT ( e(:,i), e(:,j) )
! We calculate the dipole-quadrupole and quadrupole-dipole contributions separately
! for clarity. Either part could be commented out so as to be tested separately
! together with the corresponding sections for force and torque below.
pot = 0.0
fij = 0.0
gi = 0.0
gj = 0.0
! The dipole-quadrupole potential with mu_i = 1, Q_j = 1
vij_dq = 1.5 * (ci*(1.0-5.0*cj**2)+2*cj*cij)/rij_mag**4
pot = pot + vij_dq
! The quadrupole-dipole potential with Q_i = 1, mu_j = 1
vij_qd = -1.5 * (cj*(1.0-5.0*ci**2)+2*ci*cij)/rij_mag**4
pot = pot + vij_qd
IF ( .NOT. PRESENT(f) ) RETURN
IF ( .NOT. PRESENT(t) ) THEN
WRITE ( unit=error_unit, fmt='(a)' ) 'Both f and t expected'
STOP 'Error in test_pot_dq'
END IF
IF ( ANY ( SHAPE(f) /= [3,n] ) ) THEN
WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'f shape error', SHAPE(f), 3, n
STOP 'Error in test_pot_dq'
END IF
IF ( ANY ( SHAPE(t) /= [3,n] ) ) THEN
WRITE ( unit=error_unit, fmt='(a,4i15)' ) 't shape error', SHAPE(t), 3, n
STOP 'Error in test_pot_dq'
END IF
! Forces and torques for dipole-quadrupole potential with mu_i = 1, Q_j = 1
dvdrij = -4.0*vij_dq/rij_mag
dvdci = 1.5 * (1-5.0*cj**2) / rij_mag**4
dvdcj = 3.0 * (cij-5.0*ci*cj) / rij_mag**4
dvdcij = 3.0 * cj / rij_mag**4
! Forces
fij = fij - dvdrij*sij - dvdci*(e(:,i)-ci*sij)/rij_mag - dvdcj*(e(:,j)-cj*sij)/rij_mag
! Torque gradients
gi = gi + dvdci*sij + dvdcij*e(:,j)
gj = gj + dvdcj*sij + dvdcij*e(:,i)
! Forces and torques for quadrupole-dipole potential with Q_i = 1, mu_j = 1
dvdrij = -4.0*vij_qd/rij_mag
dvdci = -3.0 * (cij-5.0*ci*cj) / rij_mag**4
dvdcj = -1.5 * (1-5.0*ci**2) / rij_mag**4
dvdcij = -3.0 * ci / rij_mag**4
! Forces
fij = fij - dvdrij*sij - dvdci*(e(:,i)-ci*sij)/rij_mag - dvdcj*(e(:,j)-cj*sij)/rij_mag
! Torque gradients
gi = gi + dvdci*sij + dvdcij*e(:,j)
gj = gj + dvdcj*sij + dvdcij*e(:,i)
! Final forces and torques
f(:,i) = fij
f(:,j) = -fij
t(:,i) = -cross_product(e(:,i),gi)
t(:,j) = -cross_product(e(:,j),gj)
END SUBROUTINE force
END MODULE test_pot_module