ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90
Revision: 2461
Committed: Mon Nov 21 22:59:02 2005 UTC (18 years, 9 months ago) by gezelter
File size: 12562 byte(s)
Log Message:
Cutoff mixing fixes have been made.

File Contents

# User Rev Content
1 gezelter 1930 !!
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 1608 !! Calculates Long Range forces Lennard-Jones interactions.
44     !! @author Charles F. Vardeman II
45     !! @author Matthew Meineke
46 gezelter 2461 !! @version $Id: LJ.F90,v 1.19 2005-11-21 22:59:01 gezelter Exp $, $Date: 2005-11-21 22:59:01 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
47 gezelter 1608
48 gezelter 1930
49 gezelter 1608 module lj
50     use atype_module
51     use vector_class
52     use simulation
53 gezelter 1628 use status
54 gezelter 1608 #ifdef IS_MPI
55     use mpiSimulation
56     #endif
57     use force_globals
58    
59     implicit none
60     PRIVATE
61 chuckv 2355 #define __FORTRAN90
62     #include "UseTheForce/DarkSide/fInteractionMap.h"
63 gezelter 2204
64 chuckv 1624 integer, parameter :: DP = selected_real_kind(15)
65 gezelter 2204
66 gezelter 2271 logical, save :: useGeometricDistanceMixing = .false.
67     logical, save :: haveMixingMap = .false.
68    
69     real(kind=DP), save :: defaultCutoff = 0.0_DP
70     logical, save :: defaultShift = .false.
71     logical, save :: haveDefaultCutoff = .false.
72    
73    
74     type, private :: LJtype
75     integer :: atid
76 gezelter 1628 real(kind=dp) :: sigma
77     real(kind=dp) :: epsilon
78 gezelter 2271 logical :: isSoftCore = .false.
79     end type LJtype
80 gezelter 2204
81 gezelter 2271 type, private :: LJList
82 chuckv 2319 integer :: Nljtypes = 0
83 gezelter 2271 integer :: currentLJtype = 0
84     type(LJtype), pointer :: LJtypes(:) => null()
85     integer, pointer :: atidToLJtype(:) => null()
86     end type LJList
87 gezelter 2204
88 gezelter 2271 type(LJList), save :: LJMap
89 gezelter 2204
90 gezelter 1628 type :: MixParameters
91     real(kind=DP) :: sigma
92     real(kind=DP) :: epsilon
93 gezelter 2271 real(kind=dp) :: sigma6
94     real(kind=dp) :: rCut
95     real(kind=dp) :: delta
96     logical :: rCutWasSet = .false.
97     logical :: shiftedPot
98     logical :: isSoftCore = .false.
99 gezelter 1628 end type MixParameters
100 gezelter 2204
101 gezelter 1628 type(MixParameters), dimension(:,:), allocatable :: MixingMap
102 gezelter 2204
103 gezelter 2271 public :: newLJtype
104     public :: setLJDefaultCutoff
105     public :: getSigma
106     public :: getEpsilon
107 gezelter 1628 public :: useGeometricMixing
108     public :: do_lj_pair
109 gezelter 2271 public :: destroyLJtypes
110 gezelter 2204
111 gezelter 1608 contains
112    
113 gezelter 2271 subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
114 gezelter 1930 integer,intent(in) :: c_ident
115 gezelter 1628 real(kind=dp),intent(in) :: sigma
116     real(kind=dp),intent(in) :: epsilon
117 gezelter 2271 integer, intent(in) :: isSoftCore
118 chuckv 1624 integer,intent(out) :: status
119 gezelter 2271 integer :: nLJTypes, ntypes, myATID
120 gezelter 1930 integer, pointer :: MatchList(:) => null()
121 gezelter 2271 integer :: current
122 gezelter 1628
123 chuckv 1624 status = 0
124 gezelter 2271 ! check to see if this is the first time into this routine...
125     if (.not.associated(LJMap%LJtypes)) then
126 gezelter 2204
127 gezelter 2271 call getMatchingElementList(atypes, "is_LennardJones", .true., &
128     nLJTypes, MatchList)
129    
130     LJMap%nLJtypes = nLJTypes
131    
132     allocate(LJMap%LJtypes(nLJTypes))
133    
134     ntypes = getSize(atypes)
135    
136     allocate(LJMap%atidToLJtype(ntypes))
137     end if
138    
139     LJMap%currentLJtype = LJMap%currentLJtype + 1
140     current = LJMap%currentLJtype
141    
142 gezelter 1930 myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
143 gezelter 2271 LJMap%atidToLJtype(myATID) = current
144     LJMap%LJtypes(current)%atid = myATID
145     LJMap%LJtypes(current)%sigma = sigma
146     LJMap%LJtypes(current)%epsilon = epsilon
147     if (isSoftCore .eq. 1) then
148     LJMap%LJtypes(current)%isSoftCore = .true.
149     else
150     LJMap%LJtypes(current)%isSoftCore = .false.
151     endif
152     end subroutine newLJtype
153 chuckv 1624
154 gezelter 2271 subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
155     real(kind=dp), intent(in) :: thisRcut
156     logical, intent(in) :: shiftedPot
157     defaultCutoff = thisRcut
158     defaultShift = shiftedPot
159     haveDefaultCutoff = .true.
160 chuckv 2319 !we only want to build LJ Mixing map if LJ is being used.
161     if(LJMap%nLJTypes /= 0) then
162     call createMixingMap()
163     end if
164 gezelter 2271 end subroutine setLJDefaultCutoff
165 gezelter 2204
166 gezelter 1633 function getSigma(atid) result (s)
167     integer, intent(in) :: atid
168 gezelter 2271 integer :: ljt1
169 gezelter 1633 real(kind=dp) :: s
170 gezelter 2204
171 gezelter 2271 if (LJMap%currentLJtype == 0) then
172     call handleError("LJ", "No members in LJMap")
173 gezelter 1633 return
174     end if
175 gezelter 2204
176 gezelter 2271 ljt1 = LJMap%atidToLJtype(atid)
177     s = LJMap%LJtypes(ljt1)%sigma
178    
179 gezelter 1633 end function getSigma
180    
181     function getEpsilon(atid) result (e)
182     integer, intent(in) :: atid
183 gezelter 2271 integer :: ljt1
184 gezelter 1633 real(kind=dp) :: e
185 gezelter 2204
186 gezelter 2271 if (LJMap%currentLJtype == 0) then
187     call handleError("LJ", "No members in LJMap")
188 gezelter 1633 return
189     end if
190 gezelter 2204
191 gezelter 2271 ljt1 = LJMap%atidToLJtype(atid)
192     e = LJMap%LJtypes(ljt1)%epsilon
193    
194 gezelter 1633 end function getEpsilon
195    
196 gezelter 1628 subroutine useGeometricMixing()
197     useGeometricDistanceMixing = .true.
198     haveMixingMap = .false.
199     return
200     end subroutine useGeometricMixing
201 gezelter 2204
202 gezelter 2271 subroutine createMixingMap()
203     integer :: nLJtypes, i, j
204     real ( kind = dp ) :: s1, s2, e1, e2
205     real ( kind = dp ) :: rcut6, tp6, tp12
206 gezelter 2289 logical :: isSoftCore1, isSoftCore2, doShift
207 gezelter 1628
208 gezelter 2271 if (LJMap%currentLJtype == 0) then
209     call handleError("LJ", "No members in LJMap")
210 gezelter 2086 return
211 gezelter 2271 end if
212 gezelter 2204
213 gezelter 2271 nLJtypes = LJMap%nLJtypes
214 gezelter 2204
215 gezelter 2086 if (.not. allocated(MixingMap)) then
216 gezelter 2271 allocate(MixingMap(nLJtypes, nLJtypes))
217 gezelter 2086 endif
218    
219 gezelter 2271 do i = 1, nLJtypes
220 gezelter 2204
221 gezelter 2271 s1 = LJMap%LJtypes(i)%sigma
222     e1 = LJMap%LJtypes(i)%epsilon
223     isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
224 gezelter 2204
225 gezelter 2271 do j = i, nLJtypes
226    
227     s2 = LJMap%LJtypes(j)%sigma
228     e2 = LJMap%LJtypes(j)%epsilon
229     isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
230    
231     MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
232 gezelter 2204
233 gezelter 2271 ! only the distance parameter uses different mixing policies
234     if (useGeometricDistanceMixing) then
235     MixingMap(i,j)%sigma = dsqrt(s1 * s2)
236     else
237     MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
238     endif
239    
240     MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
241 gezelter 2204
242 gezelter 2271 MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
243 gezelter 2204
244 gezelter 2461 if (haveDefaultCutoff) then
245     rcut6 = defaultCutoff**6
246     tp6 = MixingMap(i,j)%sigma6/rcut6
247     tp12 = tp6**2
248     MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
249     MixingMap(i,j)%shiftedPot = defaultShift
250 gezelter 2271 else
251 gezelter 2461 MixingMap(i,j)%delta = 0.0_DP
252     MixingMap(i,j)%shiftedPot = defaultShift
253     endif
254 gezelter 2204
255 gezelter 2289 if (i.ne.j) then
256     MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
257     MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
258     MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
259     MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
260     MixingMap(j,i)%delta = MixingMap(i,j)%delta
261     MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
262     MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
263     MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
264     endif
265    
266 gezelter 2271 enddo
267     enddo
268    
269 chrisfen 1754 haveMixingMap = .true.
270 gezelter 2271
271 gezelter 1628 end subroutine createMixingMap
272 gezelter 2271
273 gezelter 2461 subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
274 gezelter 1608 pot, f, do_pot)
275 gezelter 2461
276 gezelter 1608 integer, intent(in) :: atom1, atom2
277 gezelter 2271 integer :: atid1, atid2, ljt1, ljt2
278 gezelter 2461 real( kind = dp ), intent(in) :: rij, r2, rcut
279 gezelter 1608 real( kind = dp ) :: pot, sw, vpair
280     real( kind = dp ), dimension(3,nLocal) :: f
281     real( kind = dp ), intent(in), dimension(3) :: d
282     real( kind = dp ), intent(inout), dimension(3) :: fpair
283     logical, intent(in) :: do_pot
284    
285     ! local Variables
286     real( kind = dp ) :: drdx, drdy, drdz
287     real( kind = dp ) :: fx, fy, fz
288     real( kind = dp ) :: pot_temp, dudr
289     real( kind = dp ) :: sigma6
290     real( kind = dp ) :: epsilon
291 gezelter 2461 real( kind = dp ) :: r6, rc6
292     real( kind = dp ) :: t6, tp6
293     real( kind = dp ) :: t12, tp12
294 gezelter 1608 real( kind = dp ) :: delta
295 gezelter 2271 logical :: isSoftCore, shiftedPot
296 gezelter 1628 integer :: id1, id2, localError
297 gezelter 1608
298 gezelter 1628 if (.not.haveMixingMap) then
299 gezelter 2271 call createMixingMap()
300 gezelter 1628 endif
301    
302 gezelter 1608 ! Look up the correct parameters in the mixing matrix
303     #ifdef IS_MPI
304 gezelter 2271 atid1 = atid_Row(atom1)
305     atid2 = atid_Col(atom2)
306 gezelter 1608 #else
307 gezelter 2271 atid1 = atid(atom1)
308     atid2 = atid(atom2)
309 gezelter 1608 #endif
310    
311 gezelter 2271 ljt1 = LJMap%atidToLJtype(atid1)
312     ljt2 = LJMap%atidToLJtype(atid2)
313    
314     sigma6 = MixingMap(ljt1,ljt2)%sigma6
315     epsilon = MixingMap(ljt1,ljt2)%epsilon
316     delta = MixingMap(ljt1,ljt2)%delta
317     isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
318     shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
319    
320 gezelter 1608 r6 = r2 * r2 * r2
321 gezelter 2204
322 gezelter 1608 t6 = sigma6/ r6
323     t12 = t6 * t6
324 tim 2062
325 gezelter 2271 if (isSoftCore) then
326    
327     pot_temp = 4.0E0_DP * epsilon * t6
328     if (shiftedPot) then
329 gezelter 2461 rc6 = rcut**6
330     tp6 = sigma6 / rc6
331     delta =-4.0_DP*epsilon*(tp6)
332 gezelter 2271 pot_temp = pot_temp + delta
333     endif
334    
335     vpair = vpair + pot_temp
336    
337     dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
338    
339 tim 2062 else
340 gezelter 2204 pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
341 gezelter 2271 if (shiftedPot) then
342 gezelter 2461 rc6 = rcut**6
343     tp6 = sigma6 / rc6
344     tp12 = tp6*tp6
345     delta =-4.0_DP*epsilon*(tp12 - tp6)
346 gezelter 2204 pot_temp = pot_temp + delta
347     endif
348 gezelter 2271
349 gezelter 2204 vpair = vpair + pot_temp
350 gezelter 2271
351 gezelter 2204 dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
352 gezelter 1608 endif
353 gezelter 2204
354 gezelter 1608 drdx = d(1) / rij
355     drdy = d(2) / rij
356     drdz = d(3) / rij
357 gezelter 2204
358 gezelter 1608 fx = dudr * drdx
359     fy = dudr * drdy
360     fz = dudr * drdz
361 gezelter 2204
362 gezelter 1608 #ifdef IS_MPI
363     if (do_pot) then
364 gezelter 2361 pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
365     pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
366 gezelter 1608 endif
367 gezelter 2204
368 gezelter 1608 f_Row(1,atom1) = f_Row(1,atom1) + fx
369     f_Row(2,atom1) = f_Row(2,atom1) + fy
370     f_Row(3,atom1) = f_Row(3,atom1) + fz
371 gezelter 2204
372 gezelter 1608 f_Col(1,atom2) = f_Col(1,atom2) - fx
373     f_Col(2,atom2) = f_Col(2,atom2) - fy
374     f_Col(3,atom2) = f_Col(3,atom2) - fz
375 gezelter 2204
376 gezelter 1608 #else
377     if (do_pot) pot = pot + sw*pot_temp
378    
379     f(1,atom1) = f(1,atom1) + fx
380     f(2,atom1) = f(2,atom1) + fy
381     f(3,atom1) = f(3,atom1) + fz
382 gezelter 2204
383 gezelter 1608 f(1,atom2) = f(1,atom2) - fx
384     f(2,atom2) = f(2,atom2) - fy
385     f(3,atom2) = f(3,atom2) - fz
386     #endif
387 gezelter 2204
388 gezelter 1608 #ifdef IS_MPI
389     id1 = AtomRowToGlobal(atom1)
390     id2 = AtomColToGlobal(atom2)
391     #else
392     id1 = atom1
393     id2 = atom2
394     #endif
395    
396     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
397 gezelter 2204
398 gezelter 1608 fpair(1) = fpair(1) + fx
399     fpair(2) = fpair(2) + fy
400     fpair(3) = fpair(3) + fz
401    
402     endif
403    
404     return
405 gezelter 2204
406 gezelter 1608 end subroutine do_lj_pair
407 gezelter 2204
408 chuckv 2189 subroutine destroyLJTypes()
409 gezelter 2271
410     LJMap%nLJtypes = 0
411     LJMap%currentLJtype = 0
412    
413     if (associated(LJMap%LJtypes)) then
414     deallocate(LJMap%LJtypes)
415     LJMap%LJtypes => null()
416     end if
417    
418     if (associated(LJMap%atidToLJtype)) then
419     deallocate(LJMap%atidToLJtype)
420     LJMap%atidToLJtype => null()
421     end if
422    
423 chuckv 2189 haveMixingMap = .false.
424     end subroutine destroyLJTypes
425    
426 gezelter 1608 end module lj