ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/LJ.F90
Revision: 798
Committed: Wed Dec 7 19:58:18 2005 UTC (19 years, 7 months ago) by chuckv
File size: 12461 byte(s)
Log Message:
Removed geometric distance mixing from individual modules and now use
force options to detemine what the deal is.

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