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

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 gezelter 762 !! @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 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 gezelter 115 #ifdef IS_MPI
55     use mpiSimulation
56     #endif
57     use force_globals
58    
59     implicit none
60     PRIVATE
61 chuckv 656 #define __FORTRAN90
62     #include "UseTheForce/DarkSide/fInteractionMap.h"
63 gezelter 507
64 chuckv 131 integer, parameter :: DP = selected_real_kind(15)
65 gezelter 507
66 gezelter 572 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 135 real(kind=dp) :: sigma
77     real(kind=dp) :: epsilon
78 gezelter 572 logical :: isSoftCore = .false.
79     end type LJtype
80 gezelter 507
81 gezelter 572 type, private :: LJList
82 chuckv 620 integer :: Nljtypes = 0
83 gezelter 572 integer :: currentLJtype = 0
84     type(LJtype), pointer :: LJtypes(:) => null()
85     integer, pointer :: atidToLJtype(:) => null()
86     end type LJList
87 gezelter 507
88 gezelter 572 type(LJList), save :: LJMap
89 gezelter 507
90 gezelter 135 type :: MixParameters
91     real(kind=DP) :: sigma
92     real(kind=DP) :: epsilon
93 gezelter 572 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 135 end type MixParameters
100 gezelter 507
101 gezelter 135 type(MixParameters), dimension(:,:), allocatable :: MixingMap
102 gezelter 507
103 gezelter 572 public :: newLJtype
104     public :: setLJDefaultCutoff
105     public :: getSigma
106     public :: getEpsilon
107 gezelter 135 public :: useGeometricMixing
108     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 135 subroutine useGeometricMixing()
197     useGeometricDistanceMixing = .true.
198     haveMixingMap = .false.
199     return
200     end subroutine useGeometricMixing
201 gezelter 507
202 gezelter 572 subroutine createMixingMap()
203     integer :: nLJtypes, i, j
204     real ( kind = dp ) :: s1, s2, e1, e2
205     real ( kind = dp ) :: rcut6, tp6, tp12
206 gezelter 590 logical :: isSoftCore1, isSoftCore2, doShift
207 gezelter 135
208 gezelter 572 if (LJMap%currentLJtype == 0) then
209     call handleError("LJ", "No members in LJMap")
210 gezelter 402 return
211 gezelter 572 end if
212 gezelter 507
213 gezelter 572 nLJtypes = LJMap%nLJtypes
214 gezelter 507
215 gezelter 402 if (.not. allocated(MixingMap)) then
216 gezelter 572 allocate(MixingMap(nLJtypes, nLJtypes))
217 gezelter 402 endif
218    
219 gezelter 572 do i = 1, nLJtypes
220 gezelter 507
221 gezelter 572 s1 = LJMap%LJtypes(i)%sigma
222     e1 = LJMap%LJtypes(i)%epsilon
223     isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
224 gezelter 507
225 gezelter 572 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 507
233 gezelter 572 ! 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 507
242 gezelter 572 MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
243 gezelter 507
244 gezelter 762 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 572 else
251 gezelter 762 MixingMap(i,j)%delta = 0.0_DP
252     MixingMap(i,j)%shiftedPot = defaultShift
253     endif
254 gezelter 507
255 gezelter 590 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 572 enddo
267     enddo
268    
269 chrisfen 217 haveMixingMap = .true.
270 gezelter 572
271 gezelter 135 end subroutine createMixingMap
272 gezelter 572
273 gezelter 762 subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
274 gezelter 115 pot, f, do_pot)
275 gezelter 762
276 gezelter 115 integer, intent(in) :: atom1, atom2
277 gezelter 572 integer :: atid1, atid2, ljt1, ljt2
278 gezelter 762 real( kind = dp ), intent(in) :: rij, r2, rcut
279 gezelter 115 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 762 real( kind = dp ) :: r6, rc6
292     real( kind = dp ) :: t6, tp6
293     real( kind = dp ) :: t12, tp12
294 gezelter 115 real( kind = dp ) :: delta
295 gezelter 572 logical :: isSoftCore, shiftedPot
296 gezelter 135 integer :: id1, id2, localError
297 gezelter 115
298 gezelter 135 if (.not.haveMixingMap) then
299 gezelter 572 call createMixingMap()
300 gezelter 135 endif
301    
302 gezelter 115 ! Look up the correct parameters in the mixing matrix
303     #ifdef IS_MPI
304 gezelter 572 atid1 = atid_Row(atom1)
305     atid2 = atid_Col(atom2)
306 gezelter 115 #else
307 gezelter 572 atid1 = atid(atom1)
308     atid2 = atid(atom2)
309 gezelter 115 #endif
310    
311 gezelter 572 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 115 r6 = r2 * r2 * r2
321 gezelter 507
322 gezelter 115 t6 = sigma6/ r6
323     t12 = t6 * t6
324 tim 378
325 gezelter 572 if (isSoftCore) then
326    
327     pot_temp = 4.0E0_DP * epsilon * t6
328     if (shiftedPot) then
329 gezelter 762 rc6 = rcut**6
330     tp6 = sigma6 / rc6
331     delta =-4.0_DP*epsilon*(tp6)
332 gezelter 572 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 378 else
340 gezelter 507 pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
341 gezelter 572 if (shiftedPot) then
342 gezelter 762 rc6 = rcut**6
343     tp6 = sigma6 / rc6
344     tp12 = tp6*tp6
345     delta =-4.0_DP*epsilon*(tp12 - tp6)
346 gezelter 507 pot_temp = pot_temp + delta
347     endif
348 gezelter 572
349 gezelter 507 vpair = vpair + pot_temp
350 gezelter 572
351 gezelter 507 dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
352 gezelter 115 endif
353 gezelter 507
354 gezelter 115 drdx = d(1) / rij
355     drdy = d(2) / rij
356     drdz = d(3) / rij
357 gezelter 507
358 gezelter 115 fx = dudr * drdx
359     fy = dudr * drdy
360     fz = dudr * drdz
361 gezelter 507
362 gezelter 115 #ifdef IS_MPI
363     if (do_pot) then
364 gezelter 662 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 115 endif
367 gezelter 507
368 gezelter 115 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 507
372 gezelter 115 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 507
376 gezelter 115 #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 507
383 gezelter 115 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 507
388 gezelter 115 #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 507
398 gezelter 115 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 507
406 gezelter 115 end subroutine do_lj_pair
407 gezelter 507
408 chuckv 492 subroutine destroyLJTypes()
409 gezelter 572
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 492 haveMixingMap = .false.
424     end subroutine destroyLJTypes
425    
426 gezelter 115 end module lj