ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/LJ.F90
Revision: 537
Committed: Thu May 19 15:49:53 2005 UTC (20 years, 2 months ago) by tim
File size: 12926 byte(s)
Log Message:
fix bug in NPAT and NPrT

File Contents

# User Rev Content
1 gezelter 246 !!
2     !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     !!
4     !! The University of Notre Dame grants you ("Licensee") a
5     !! non-exclusive, royalty free, license to use, modify and
6     !! redistribute this software in source and binary code form, provided
7     !! that the following conditions are met:
8     !!
9     !! 1. Acknowledgement of the program authors must be made in any
10     !! publication of scientific results based in part on use of the
11     !! program. An acceptable form of acknowledgement is citation of
12     !! the article in which the program was described (Matthew
13     !! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     !! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     !! Parallel Simulation Engine for Molecular Dynamics,"
16     !! J. Comput. Chem. 26, pp. 252-271 (2005))
17     !!
18     !! 2. Redistributions of source code must retain the above copyright
19     !! notice, this list of conditions and the following disclaimer.
20     !!
21     !! 3. Redistributions in binary form must reproduce the above copyright
22     !! notice, this list of conditions and the following disclaimer in the
23     !! documentation and/or other materials provided with the
24     !! distribution.
25     !!
26     !! This software is provided "AS IS," without a warranty of any
27     !! kind. All express or implied conditions, representations and
28     !! warranties, including any implied warranty of merchantability,
29     !! fitness for a particular purpose or non-infringement, are hereby
30     !! excluded. The University of Notre Dame and its licensors shall not
31     !! be liable for any damages suffered by licensee as a result of
32     !! using, modifying or distributing the software or its
33     !! derivatives. In no event will the University of Notre Dame or its
34     !! licensors be liable for any lost revenue, profit or data, or for
35     !! direct, indirect, special, consequential, incidental or punitive
36     !! damages, however caused and regardless of the theory of liability,
37     !! arising out of the use of or inability to use software, even if the
38     !! University of Notre Dame has been advised of the possibility of
39     !! such damages.
40     !!
41    
42    
43 gezelter 115 !! Calculates Long Range forces Lennard-Jones interactions.
44     !! @author Charles F. Vardeman II
45     !! @author Matthew Meineke
46 tim 537 !! @version $Id: LJ.F90,v 1.12 2005-05-19 15:49:53 tim Exp $, $Date: 2005-05-19 15:49:53 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $
47 gezelter 115
48 gezelter 246
49 gezelter 115 module lj
50     use atype_module
51     use switcheroo
52     use vector_class
53     use simulation
54 gezelter 135 use status
55 gezelter 115 #ifdef IS_MPI
56     use mpiSimulation
57     #endif
58     use force_globals
59    
60     implicit none
61     PRIVATE
62 gezelter 507
63 chuckv 131 integer, parameter :: DP = selected_real_kind(15)
64 gezelter 507
65 gezelter 135 type, private :: LjType
66 gezelter 246 integer :: c_ident
67     integer :: atid
68 gezelter 135 real(kind=dp) :: sigma
69     real(kind=dp) :: epsilon
70 tim 378 logical :: soft_pot
71 gezelter 135 end type LjType
72 gezelter 507
73 gezelter 135 type(LjType), dimension(:), allocatable :: ParameterMap
74 gezelter 507
75 gezelter 135 logical, save :: haveMixingMap = .false.
76 gezelter 507
77 gezelter 135 type :: MixParameters
78     real(kind=DP) :: sigma
79     real(kind=DP) :: epsilon
80     real(kind=dp) :: sigma6
81     real(kind=dp) :: tp6
82     real(kind=dp) :: tp12
83     real(kind=dp) :: delta
84 tim 378 logical :: soft_pot
85 gezelter 135 end type MixParameters
86 gezelter 507
87 gezelter 135 type(MixParameters), dimension(:,:), allocatable :: MixingMap
88 gezelter 507
89 gezelter 135 real(kind=DP), save :: LJ_rcut
90     logical, save :: have_rcut = .false.
91     logical, save :: LJ_do_shift = .false.
92     logical, save :: useGeometricDistanceMixing = .false.
93 gezelter 507
94 gezelter 135 !! Public methods and data
95 gezelter 507
96 gezelter 135 public :: setCutoffLJ
97     public :: useGeometricMixing
98     public :: do_lj_pair
99     public :: newLJtype
100 gezelter 140 public :: getSigma
101     public :: getEpsilon
102 chuckv 492 public :: destroyLJTypes
103 gezelter 507
104 gezelter 115 contains
105    
106 tim 378 subroutine newLJtype(c_ident, sigma, epsilon, soft_pot, status)
107 gezelter 246 integer,intent(in) :: c_ident
108 gezelter 135 real(kind=dp),intent(in) :: sigma
109     real(kind=dp),intent(in) :: epsilon
110 tim 378 integer, intent(in) :: soft_pot
111 chuckv 131 integer,intent(out) :: status
112 gezelter 246 integer :: nATypes, myATID
113     integer, pointer :: MatchList(:) => null()
114 gezelter 135
115 chuckv 131 status = 0
116 gezelter 507
117 gezelter 246 myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
118 chuckv 131
119 gezelter 135 if (.not.allocated(ParameterMap)) then
120 gezelter 507
121 gezelter 246 !call getMatchingElementList(atypes, "is_LennardJones", .true., &
122     ! nLJTypes, MatchList)
123 gezelter 135 nAtypes = getSize(atypes)
124 chuckv 131 if (nAtypes == 0) then
125 gezelter 135 status = -1
126     return
127 chuckv 131 end if
128 gezelter 507
129 gezelter 135 if (.not. allocated(ParameterMap)) then
130     allocate(ParameterMap(nAtypes))
131     endif
132 gezelter 507
133 chuckv 131 end if
134    
135 gezelter 246 if (myATID .gt. size(ParameterMap)) then
136 gezelter 135 status = -1
137     return
138     endif
139 gezelter 507
140 gezelter 135 ! set the values for ParameterMap for this atom type:
141    
142 gezelter 246 ParameterMap(myATID)%c_ident = c_ident
143     ParameterMap(myATID)%atid = myATID
144     ParameterMap(myATID)%epsilon = epsilon
145     ParameterMap(myATID)%sigma = sigma
146 tim 378 if (soft_pot == 1) then
147 gezelter 507 ParameterMap(myATID)%soft_pot = .true.
148 tim 378 else
149 gezelter 507 ParameterMap(myATID)%soft_pot = .false.
150 tim 378 endif
151 chuckv 131 end subroutine newLJtype
152 gezelter 140
153     function getSigma(atid) result (s)
154     integer, intent(in) :: atid
155     integer :: localError
156     real(kind=dp) :: s
157 gezelter 507
158 gezelter 140 if (.not.allocated(ParameterMap)) then
159     call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
160     return
161     end if
162 gezelter 507
163 gezelter 140 s = ParameterMap(atid)%sigma
164     end function getSigma
165    
166     function getEpsilon(atid) result (e)
167     integer, intent(in) :: atid
168     integer :: localError
169     real(kind=dp) :: e
170 gezelter 507
171 gezelter 140 if (.not.allocated(ParameterMap)) then
172 gezelter 246 call handleError("LJ", "no ParameterMap was present before first call of getEpsilon!")
173 gezelter 140 return
174     end if
175 gezelter 507
176 gezelter 140 e = ParameterMap(atid)%epsilon
177     end function getEpsilon
178    
179    
180 gezelter 115 subroutine setCutoffLJ(rcut, do_shift, status)
181     logical, intent(in):: do_shift
182     integer :: status, myStatus
183     real(kind=dp) :: rcut
184    
185     #define __FORTRAN90
186     #include "UseTheForce/fSwitchingFunction.h"
187    
188     status = 0
189    
190     LJ_rcut = rcut
191     LJ_do_shift = do_shift
192     call set_switch(LJ_SWITCH, rcut, rcut)
193 gezelter 135 have_rcut = .true.
194 gezelter 507
195 gezelter 115 return
196     end subroutine setCutoffLJ
197 gezelter 135
198     subroutine useGeometricMixing()
199     useGeometricDistanceMixing = .true.
200     haveMixingMap = .false.
201     return
202     end subroutine useGeometricMixing
203 gezelter 507
204 gezelter 135 subroutine createMixingMap(status)
205 gezelter 246 integer :: nATIDs
206 gezelter 115 integer :: status
207     integer :: i
208     integer :: j
209 gezelter 135 real ( kind = dp ) :: Sigma_i, Sigma_j
210     real ( kind = dp ) :: Epsilon_i, Epsilon_j
211 gezelter 115 real ( kind = dp ) :: rcut6
212 gezelter 402 logical :: i_is_LJ, j_is_LJ
213 gezelter 135
214 gezelter 115 status = 0
215 gezelter 402
216     if (.not. allocated(ParameterMap)) then
217     call handleError("LJ", "no ParameterMap was present before call of createMixingMap!")
218     status = -1
219     return
220     endif
221 gezelter 507
222 gezelter 246 nATIDs = size(ParameterMap)
223 gezelter 507
224 gezelter 246 if (nATIDs == 0) then
225 gezelter 115 status = -1
226     return
227     end if
228 gezelter 135
229 gezelter 402 if (.not. allocated(MixingMap)) then
230     allocate(MixingMap(nATIDs, nATIDs))
231     endif
232    
233 gezelter 135 if (.not.have_rcut) then
234     status = -1
235     return
236 gezelter 115 endif
237 gezelter 507
238 gezelter 115 rcut6 = LJ_rcut**6
239 gezelter 507
240 gezelter 135 ! This loops through all atypes, even those that don't support LJ forces.
241 gezelter 246 do i = 1, nATIDs
242 gezelter 402 call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
243     if (i_is_LJ) then
244     Epsilon_i = ParameterMap(i)%epsilon
245     Sigma_i = ParameterMap(i)%sigma
246 gezelter 507
247 gezelter 402 ! do self mixing rule
248     MixingMap(i,i)%sigma = Sigma_i
249     MixingMap(i,i)%sigma6 = Sigma_i ** 6
250     MixingMap(i,i)%tp6 = (MixingMap(i,i)%sigma6)/rcut6
251     MixingMap(i,i)%tp12 = (MixingMap(i,i)%tp6) ** 2
252     MixingMap(i,i)%epsilon = Epsilon_i
253     MixingMap(i,i)%delta = -4.0_DP * MixingMap(i,i)%epsilon * &
254     (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
255     MixingMap(i,i)%soft_pot = ParameterMap(i)%soft_pot
256 gezelter 507
257 gezelter 402 do j = i + 1, nATIDs
258     call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
259 gezelter 507
260 gezelter 402 if (j_is_LJ) then
261     Epsilon_j = ParameterMap(j)%epsilon
262     Sigma_j = ParameterMap(j)%sigma
263 gezelter 507
264 gezelter 402 ! only the distance parameter uses different mixing policies
265     if (useGeometricDistanceMixing) then
266     ! only for OPLS as far as we can tell
267     MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
268     else
269     ! everyone else
270     MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
271     endif
272 gezelter 507
273 gezelter 402 ! energy parameter is always geometric mean:
274     MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
275 gezelter 507
276 gezelter 402 MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
277     MixingMap(i,j)%tp6 = MixingMap(i,j)%sigma6/rcut6
278     MixingMap(i,j)%tp12 = (MixingMap(i,j)%tp6) ** 2
279 gezelter 507
280 gezelter 402 MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
281     (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
282 gezelter 507
283 gezelter 402 MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
284 gezelter 507
285    
286 gezelter 402 MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
287     MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
288     MixingMap(j,i)%tp6 = MixingMap(i,j)%tp6
289     MixingMap(j,i)%tp12 = MixingMap(i,j)%tp12
290     MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
291     MixingMap(j,i)%delta = MixingMap(i,j)%delta
292     MixingMap(j,i)%soft_pot = MixingMap(i,j)%soft_pot
293     endif
294     end do
295     endif
296 gezelter 115 end do
297 gezelter 507
298 chrisfen 217 haveMixingMap = .true.
299    
300 gezelter 135 end subroutine createMixingMap
301 gezelter 507
302 gezelter 115 subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
303     pot, f, do_pot)
304    
305     integer, intent(in) :: atom1, atom2
306     real( kind = dp ), intent(in) :: rij, r2
307     real( kind = dp ) :: pot, sw, vpair
308     real( kind = dp ), dimension(3,nLocal) :: f
309     real( kind = dp ), intent(in), dimension(3) :: d
310     real( kind = dp ), intent(inout), dimension(3) :: fpair
311     logical, intent(in) :: do_pot
312    
313     ! local Variables
314     real( kind = dp ) :: drdx, drdy, drdz
315     real( kind = dp ) :: fx, fy, fz
316     real( kind = dp ) :: pot_temp, dudr
317     real( kind = dp ) :: sigma6
318     real( kind = dp ) :: epsilon
319     real( kind = dp ) :: r6
320     real( kind = dp ) :: t6
321     real( kind = dp ) :: t12
322     real( kind = dp ) :: delta
323 tim 378 logical :: soft_pot
324 gezelter 135 integer :: id1, id2, localError
325 gezelter 115
326 gezelter 135 if (.not.haveMixingMap) then
327     localError = 0
328     call createMixingMap(localError)
329     if ( localError .ne. 0 ) then
330     call handleError("LJ", "MixingMap creation failed!")
331     return
332     end if
333     endif
334    
335 gezelter 115 ! Look up the correct parameters in the mixing matrix
336     #ifdef IS_MPI
337 gezelter 135 sigma6 = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
338     epsilon = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
339     delta = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
340 tim 378 soft_pot = MixingMap(atid_Row(atom1),atid_Col(atom2))%soft_pot
341 gezelter 115 #else
342 gezelter 135 sigma6 = MixingMap(atid(atom1),atid(atom2))%sigma6
343     epsilon = MixingMap(atid(atom1),atid(atom2))%epsilon
344     delta = MixingMap(atid(atom1),atid(atom2))%delta
345 tim 378 soft_pot = MixingMap(atid(atom1),atid(atom2))%soft_pot
346 gezelter 115 #endif
347    
348     r6 = r2 * r2 * r2
349 gezelter 507
350 gezelter 115 t6 = sigma6/ r6
351     t12 = t6 * t6
352 tim 378
353     if (soft_pot) then
354    
355 tim 537 pot_temp = 4.0E0_DP * epsilon * t6
356     if (LJ_do_shift) then
357     pot_temp = pot_temp + delta
358     endif
359 gezelter 507
360 tim 537 vpair = vpair + pot_temp
361 gezelter 507
362 tim 537 dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
363    
364 tim 378 else
365 gezelter 507 pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
366     if (LJ_do_shift) then
367     pot_temp = pot_temp + delta
368     endif
369    
370     vpair = vpair + pot_temp
371    
372     dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
373 gezelter 115 endif
374 gezelter 507
375 gezelter 115 drdx = d(1) / rij
376     drdy = d(2) / rij
377     drdz = d(3) / rij
378 gezelter 507
379 gezelter 115 fx = dudr * drdx
380     fy = dudr * drdy
381     fz = dudr * drdz
382 gezelter 507
383    
384 gezelter 115 #ifdef IS_MPI
385     if (do_pot) then
386     pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
387     pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
388     endif
389 gezelter 507
390 gezelter 115 f_Row(1,atom1) = f_Row(1,atom1) + fx
391     f_Row(2,atom1) = f_Row(2,atom1) + fy
392     f_Row(3,atom1) = f_Row(3,atom1) + fz
393 gezelter 507
394 gezelter 115 f_Col(1,atom2) = f_Col(1,atom2) - fx
395     f_Col(2,atom2) = f_Col(2,atom2) - fy
396     f_Col(3,atom2) = f_Col(3,atom2) - fz
397 gezelter 507
398 gezelter 115 #else
399     if (do_pot) pot = pot + sw*pot_temp
400    
401     f(1,atom1) = f(1,atom1) + fx
402     f(2,atom1) = f(2,atom1) + fy
403     f(3,atom1) = f(3,atom1) + fz
404 gezelter 507
405 gezelter 115 f(1,atom2) = f(1,atom2) - fx
406     f(2,atom2) = f(2,atom2) - fy
407     f(3,atom2) = f(3,atom2) - fz
408     #endif
409 gezelter 507
410 gezelter 115 #ifdef IS_MPI
411     id1 = AtomRowToGlobal(atom1)
412     id2 = AtomColToGlobal(atom2)
413     #else
414     id1 = atom1
415     id2 = atom2
416     #endif
417    
418     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
419 gezelter 507
420 gezelter 115 fpair(1) = fpair(1) + fx
421     fpair(2) = fpair(2) + fy
422     fpair(3) = fpair(3) + fz
423    
424     endif
425    
426     return
427 gezelter 507
428 gezelter 115 end subroutine do_lj_pair
429 gezelter 507
430 chuckv 492 subroutine destroyLJTypes()
431     if(allocated(ParameterMap)) deallocate(ParameterMap)
432     if(allocated(MixingMap)) deallocate(MixingMap)
433     haveMixingMap = .false.
434     end subroutine destroyLJTypes
435    
436 gezelter 507
437 gezelter 115 end module lj