← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to tools/i-pi/drivers/LJ.f90

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/i-pi/drivers/LJ.f90
 
1
! This performs the calculations necessary to run a Lennard-Jones (LJ) 
 
2
! simulation.
 
3
! 
 
4
! Copyright (C) 2013, Joshua More and Michele Ceriotti
 
5
! 
 
6
! Permission is hereby granted, free of charge, to any person obtaining
 
7
! a copy of this software and associated documentation files (the
 
8
! "Software"), to deal in the Software without restriction, including
 
9
! without limitation the rights to use, copy, modify, merge, publish,
 
10
! distribute, sublicense, and/or sell copies of the Software, and to
 
11
! permit persons to whom the Software is furnished to do so, subject to
 
12
! the following conditions:
 
13
! 
 
14
! The above copyright notice and this permission notice shall be included
 
15
! in all copies or substantial portions of the Software.
 
16
! 
 
17
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 
18
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 
19
! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 
20
! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 
21
! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 
22
! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 
23
! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
24
!
 
25
!
 
26
! This contains the functions that calculate the potential, forces and
 
27
! virial tensor of a single-component LJ system.
 
28
! Includes functions which calculate the long-range correction terms for a
 
29
! simulation with a sharp nearest-neighbour cut-off.
 
30
!
 
31
! Functions:
 
32
!    LJ_functions: Calculates the LJ pair potential and the magnitude of the
 
33
!       forces acting on a pair of atoms.
 
34
!    LJ_fij: Calculates the LJ pair potential and force vector for the
 
35
!       interaction of a pair of atoms.
 
36
!    LJ_longrange: Calculates the long range correction to the potential
 
37
!       and virial.
 
38
!    LJ_getall: Calculates the potential and virial of the system and the
 
39
!       forces acting on all the atoms.
 
40
 
 
41
      MODULE LJ
 
42
         USE DISTANCE
 
43
      IMPLICIT NONE
 
44
 
 
45
      DOUBLE PRECISION, PARAMETER :: four_tau_by_3 = 8.3775804095727811d0 
 
46
 
 
47
      CONTAINS
 
48
 
 
49
         SUBROUTINE LJ_functions(sigma, eps, r, pot, force)
 
50
            ! Calculates the magnitude of the LJ force and potential between
 
51
            ! a pair of atoms at a given distance from each other.
 
52
            !
 
53
            ! Args:
 
54
            !    sigma: The LJ distance parameter.
 
55
            !    eps: The LJ energy parameter.
 
56
            !    r: The separation of the atoms.
 
57
            !    pot: The LJ interaction potential.
 
58
            !    force: The magnitude of the LJ force.
 
59
 
 
60
            DOUBLE PRECISION, INTENT(IN) :: sigma
 
61
            DOUBLE PRECISION, INTENT(IN) :: eps
 
62
            DOUBLE PRECISION, INTENT(IN) :: r
 
63
            DOUBLE PRECISION, INTENT(OUT) :: pot
 
64
            DOUBLE PRECISION, INTENT(OUT) :: force
 
65
 
 
66
            DOUBLE PRECISION sigma_by_r6
 
67
 
 
68
            sigma_by_r6 = sigma/r
 
69
            sigma_by_r6 = sigma_by_r6*sigma_by_r6*sigma_by_r6
 
70
            sigma_by_r6 = sigma_by_r6*sigma_by_r6
 
71
 
 
72
            pot = 4*eps*(sigma_by_r6*(sigma_by_r6 - 1))
 
73
            force = 24*eps*(sigma_by_r6*(2*sigma_by_r6 - 1)/r)
 
74
 
 
75
         END SUBROUTINE
 
76
 
 
77
         SUBROUTINE LJ_fij(sigma, eps, rij, r, pot, fij)
 
78
            ! This calculates the LJ potential energy and the magnitude and
 
79
            ! direction of the force acting on a pair of atoms.
 
80
            !
 
81
            ! Args:
 
82
            !    sigma: The LJ distance parameter.
 
83
            !    eps: The LJ energy parameter.
 
84
            !    rij: The vector joining the two atoms.
 
85
            !    r: The separation of the two atoms.
 
86
            !    pot: The LJ interaction potential.
 
87
            !    fij: The LJ force vector.
 
88
 
 
89
            DOUBLE PRECISION, INTENT(IN) :: sigma
 
90
            DOUBLE PRECISION, INTENT(IN) :: eps
 
91
            DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: rij
 
92
            DOUBLE PRECISION, INTENT(IN) :: r
 
93
            DOUBLE PRECISION, INTENT(OUT) :: pot
 
94
            DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: fij
 
95
   
 
96
            DOUBLE PRECISION f_tot
 
97
   
 
98
            CALL LJ_functions(sigma, eps, r, pot, f_tot)
 
99
            fij = f_tot*rij/r
 
100
   
 
101
         END SUBROUTINE
 
102
 
 
103
         SUBROUTINE LJ_longrange(rc, sigma, eps, natoms, volume, pot_lr, vir_lr)
 
104
            ! Calculates the long range correction to the total potential and
 
105
            ! virial pressure.
 
106
            !
 
107
            ! Uses the tail correction for a sharp cut-off, with no smoothing
 
108
            ! function, as derived in Martyna and Hughes, Journal of Chemical
 
109
            ! Physics, 110, 3275, (1999).
 
110
            !
 
111
            ! Args:
 
112
            !    rc: The cut-off radius.
 
113
            !    sigma: The LJ distance parameter.
 
114
            !    eps: The LJ energy parameter.
 
115
            !    natoms: The number of atoms in the system.
 
116
            !    volume: The volume of the system box.
 
117
            !    pot_lr: The tail correction to the LJ interaction potential.
 
118
            !    vir_lr: The tail correction to the LJ virial pressure.
 
119
 
 
120
            DOUBLE PRECISION, INTENT(IN) :: rc
 
121
            DOUBLE PRECISION, INTENT(IN) :: sigma
 
122
            DOUBLE PRECISION, INTENT(IN) :: eps
 
123
            INTEGER, INTENT(IN) :: natoms
 
124
            DOUBLE PRECISION, INTENT(IN) :: volume
 
125
            DOUBLE PRECISION, INTENT(OUT) :: pot_lr
 
126
            DOUBLE PRECISION, INTENT(OUT) :: vir_lr
 
127
 
 
128
            DOUBLE PRECISION sbyr, s3byr3, s6byr3, s6byr6, prefactor
 
129
 
 
130
            sbyr = sigma/rc
 
131
            s3byr3 = sbyr*sbyr*sbyr
 
132
            s6byr6 = s3byr3*s3byr3
 
133
            prefactor = four_tau_by_3*natoms*natoms*eps/volume
 
134
            prefactor = prefactor*s3byr3*sigma*sigma*sigma
 
135
 
 
136
            pot_lr = prefactor*(s6byr6/3-1)
 
137
            vir_lr = prefactor*(s6byr6-1) + pot_lr
 
138
 
 
139
         END SUBROUTINE
 
140
 
 
141
         SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial)
 
142
            ! Calculates the LJ potential energy and virial and the forces 
 
143
            ! acting on all the atoms.
 
144
            !
 
145
            ! Args:
 
146
            !    rc: The cut-off radius.
 
147
            !    sigma: The LJ distance parameter.
 
148
            !    eps: The LJ energy parameter.
 
149
            !    natoms: The number of atoms in the system.
 
150
            !    atoms: A vector holding all the atom positions.
 
151
            !    cell_h: The simulation box cell vector matrix.
 
152
            !    cell_ih: The inverse of the simulation box cell vector matrix.
 
153
            !    index_list: A array giving the last index of n_list that 
 
154
            !       gives the neighbours of a given atom.
 
155
            !    n_list: An array giving the indices of the atoms that neighbour
 
156
            !       the atom determined by index_list.
 
157
            !    pot: The total potential energy of the system.
 
158
            !    forces: An array giving the forces acting on all the atoms.
 
159
            !    virial: The virial tensor, not divided by the volume.
 
160
 
 
161
            DOUBLE PRECISION, INTENT(IN) :: rc
 
162
            DOUBLE PRECISION, INTENT(IN) :: sigma
 
163
            DOUBLE PRECISION, INTENT(IN) :: eps
 
164
            INTEGER, INTENT(IN) :: natoms
 
165
            DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(IN) :: atoms
 
166
            DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_h
 
167
            DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_ih
 
168
            INTEGER, DIMENSION(natoms), INTENT(IN) :: index_list
 
169
            INTEGER, DIMENSION(natoms*(natoms-1)/2), INTENT(IN) :: n_list
 
170
            DOUBLE PRECISION, INTENT(OUT) :: pot
 
171
            DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(OUT) :: forces
 
172
            DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: virial
 
173
 
 
174
            INTEGER i, j, k, l, start
 
175
            DOUBLE PRECISION, DIMENSION(3) :: fij, rij
 
176
            DOUBLE PRECISION r2, pot_ij, pot_lr, vir_lr, volume
 
177
 
 
178
            forces = 0.0d0
 
179
            pot = 0.0d0
 
180
            virial = 0.0d0
 
181
 
 
182
            start = 1
 
183
 
 
184
            DO i = 1, natoms - 1
 
185
               ! Only loops over the neigbour list, not all the atoms.
 
186
               DO j = start, index_list(i)
 
187
                  CALL vector_separation(cell_h, cell_ih, atoms(i,:), atoms(n_list(j),:), rij, r2)
 
188
                  IF (r2 < rc*rc) THEN ! Only calculates contributions between neighbouring particles.
 
189
                     CALL LJ_fij(sigma, eps, rij, sqrt(r2), pot_ij, fij)
 
190
 
 
191
                     forces(i,:) = forces(i,:) + fij
 
192
                     forces(n_list(j),:) = forces(n_list(j),:) - fij
 
193
                     pot = pot + pot_ij
 
194
                     DO k = 1, 3
 
195
                        DO l = k, 3
 
196
                           ! Only the upper triangular elements calculated.
 
197
                           virial(k,l) = virial(k,l) + fij(k)*rij(l)
 
198
                        ENDDO
 
199
                     ENDDO
 
200
                  ENDIF
 
201
               ENDDO
 
202
               start = index_list(i) + 1
 
203
            ENDDO
 
204
 
 
205
            ! Assuming an upper-triangular vector matrix for the simulation box.
 
206
            volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
 
207
            CALL LJ_longrange(rc, sigma, eps, natoms, volume, pot_lr, vir_lr)
 
208
            pot = pot + pot_lr
 
209
            DO k = 1, 3
 
210
               virial(k,k) = virial(k,k) + vir_lr
 
211
            ENDDO
 
212
 
 
213
         END SUBROUTINE
 
214
 
 
215
      END MODULE