ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/LJ.F90
Revision: 662
Committed: Wed Oct 12 21:00:50 2005 UTC (19 years, 9 months ago) by gezelter
File size: 14082 byte(s)
Log Message:
simplifying long range potential array

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 662 !! @version $Id: LJ.F90,v 1.18 2005-10-12 21:00:49 gezelter Exp $, $Date: 2005-10-12 21:00:49 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
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 :: setLJUniformCutoff
106     public :: setLJCutoffByTypes
107     public :: getSigma
108     public :: getEpsilon
109 gezelter 135 public :: useGeometricMixing
110     public :: do_lj_pair
111 gezelter 572 public :: destroyLJtypes
112 gezelter 507
113 gezelter 115 contains
114    
115 gezelter 572 subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
116 gezelter 246 integer,intent(in) :: c_ident
117 gezelter 135 real(kind=dp),intent(in) :: sigma
118     real(kind=dp),intent(in) :: epsilon
119 gezelter 572 integer, intent(in) :: isSoftCore
120 chuckv 131 integer,intent(out) :: status
121 gezelter 572 integer :: nLJTypes, ntypes, myATID
122 gezelter 246 integer, pointer :: MatchList(:) => null()
123 gezelter 572 integer :: current
124 gezelter 135
125 chuckv 131 status = 0
126 gezelter 572 ! check to see if this is the first time into this routine...
127     if (.not.associated(LJMap%LJtypes)) then
128 gezelter 507
129 gezelter 572 call getMatchingElementList(atypes, "is_LennardJones", .true., &
130     nLJTypes, MatchList)
131    
132     LJMap%nLJtypes = nLJTypes
133    
134     allocate(LJMap%LJtypes(nLJTypes))
135    
136     ntypes = getSize(atypes)
137    
138     allocate(LJMap%atidToLJtype(ntypes))
139     end if
140    
141     LJMap%currentLJtype = LJMap%currentLJtype + 1
142     current = LJMap%currentLJtype
143    
144 gezelter 246 myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
145 gezelter 572 LJMap%atidToLJtype(myATID) = current
146     LJMap%LJtypes(current)%atid = myATID
147     LJMap%LJtypes(current)%sigma = sigma
148     LJMap%LJtypes(current)%epsilon = epsilon
149     if (isSoftCore .eq. 1) then
150     LJMap%LJtypes(current)%isSoftCore = .true.
151     else
152     LJMap%LJtypes(current)%isSoftCore = .false.
153     endif
154     end subroutine newLJtype
155 chuckv 131
156 gezelter 572 subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
157     real(kind=dp), intent(in) :: thisRcut
158     logical, intent(in) :: shiftedPot
159     defaultCutoff = thisRcut
160     defaultShift = shiftedPot
161     haveDefaultCutoff = .true.
162 chuckv 620 !we only want to build LJ Mixing map if LJ is being used.
163     if(LJMap%nLJTypes /= 0) then
164     call createMixingMap()
165     end if
166 gezelter 572 end subroutine setLJDefaultCutoff
167 gezelter 507
168 gezelter 572 subroutine setLJUniformCutoff(thisRcut, shiftedPot)
169     real(kind=dp), intent(in) :: thisRcut
170     logical, intent(in) :: shiftedPot
171     integer :: nLJtypes, i, j
172 gezelter 507
173 gezelter 572 if (LJMap%currentLJtype == 0) then
174     call handleError("LJ", "No members in LJMap")
175     return
176 chuckv 131 end if
177    
178 gezelter 572 nLJtypes = LJMap%nLJtypes
179     if (.not. allocated(MixingMap)) then
180     allocate(MixingMap(nLJtypes, nLJtypes))
181 gezelter 135 endif
182 gezelter 507
183 gezelter 572 do i = 1, nLJtypes
184     do j = 1, nLJtypes
185     MixingMap(i,j)%rCut = thisRcut
186     MixingMap(i,j)%shiftedPot = shiftedPot
187     MixingMap(i,j)%rCutWasSet = .true.
188     enddo
189     enddo
190     call createMixingMap()
191     end subroutine setLJUniformCutoff
192 gezelter 135
193 gezelter 572 subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
194     integer, intent(in) :: atid1, atid2
195     real(kind=dp), intent(in) :: thisRcut
196     logical, intent(in) :: shiftedPot
197     integer :: nLJtypes, ljt1, ljt2
198    
199     if (LJMap%currentLJtype == 0) then
200     call handleError("LJ", "No members in LJMap")
201     return
202     end if
203    
204     nLJtypes = LJMap%nLJtypes
205     if (.not. allocated(MixingMap)) then
206     allocate(MixingMap(nLJtypes, nLJtypes))
207 tim 378 endif
208 gezelter 140
209 gezelter 572 ljt1 = LJMap%atidToLJtype(atid1)
210     ljt2 = LJMap%atidToLJtype(atid2)
211    
212     MixingMap(ljt1,ljt2)%rCut = thisRcut
213     MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
214     MixingMap(ljt1,ljt2)%rCutWasSet = .true.
215     MixingMap(ljt2,ljt1)%rCut = thisRcut
216     MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
217     MixingMap(ljt2,ljt1)%rCutWasSet = .true.
218    
219     call createMixingMap()
220     end subroutine setLJCutoffByTypes
221    
222 gezelter 140 function getSigma(atid) result (s)
223     integer, intent(in) :: atid
224 gezelter 572 integer :: ljt1
225 gezelter 140 real(kind=dp) :: s
226 gezelter 507
227 gezelter 572 if (LJMap%currentLJtype == 0) then
228     call handleError("LJ", "No members in LJMap")
229 gezelter 140 return
230     end if
231 gezelter 507
232 gezelter 572 ljt1 = LJMap%atidToLJtype(atid)
233     s = LJMap%LJtypes(ljt1)%sigma
234    
235 gezelter 140 end function getSigma
236    
237     function getEpsilon(atid) result (e)
238     integer, intent(in) :: atid
239 gezelter 572 integer :: ljt1
240 gezelter 140 real(kind=dp) :: e
241 gezelter 507
242 gezelter 572 if (LJMap%currentLJtype == 0) then
243     call handleError("LJ", "No members in LJMap")
244 gezelter 140 return
245     end if
246 gezelter 507
247 gezelter 572 ljt1 = LJMap%atidToLJtype(atid)
248     e = LJMap%LJtypes(ljt1)%epsilon
249    
250 gezelter 140 end function getEpsilon
251    
252 gezelter 135 subroutine useGeometricMixing()
253     useGeometricDistanceMixing = .true.
254     haveMixingMap = .false.
255     return
256     end subroutine useGeometricMixing
257 gezelter 507
258 gezelter 572 subroutine createMixingMap()
259     integer :: nLJtypes, i, j
260     real ( kind = dp ) :: s1, s2, e1, e2
261     real ( kind = dp ) :: rcut6, tp6, tp12
262 gezelter 590 logical :: isSoftCore1, isSoftCore2, doShift
263 gezelter 135
264 gezelter 572 if (LJMap%currentLJtype == 0) then
265     call handleError("LJ", "No members in LJMap")
266 gezelter 402 return
267 gezelter 572 end if
268 gezelter 507
269 gezelter 572 nLJtypes = LJMap%nLJtypes
270 gezelter 507
271 gezelter 402 if (.not. allocated(MixingMap)) then
272 gezelter 572 allocate(MixingMap(nLJtypes, nLJtypes))
273 gezelter 402 endif
274    
275 gezelter 572 do i = 1, nLJtypes
276 gezelter 507
277 gezelter 572 s1 = LJMap%LJtypes(i)%sigma
278     e1 = LJMap%LJtypes(i)%epsilon
279     isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
280 gezelter 507
281 gezelter 572 do j = i, nLJtypes
282    
283     s2 = LJMap%LJtypes(j)%sigma
284     e2 = LJMap%LJtypes(j)%epsilon
285     isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
286    
287     MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
288 gezelter 507
289 gezelter 572 ! only the distance parameter uses different mixing policies
290     if (useGeometricDistanceMixing) then
291     MixingMap(i,j)%sigma = dsqrt(s1 * s2)
292     else
293     MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
294     endif
295    
296     MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
297 gezelter 507
298 gezelter 572 MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
299 gezelter 507
300 gezelter 572 if (MixingMap(i,j)%rCutWasSet) then
301     rcut6 = (MixingMap(i,j)%rcut)**6
302     else
303     if (haveDefaultCutoff) then
304     rcut6 = defaultCutoff**6
305 gezelter 590 doShift = defaultShift
306 gezelter 572 else
307     call handleError("LJ", "No specified or default cutoff value!")
308 gezelter 402 endif
309 gezelter 572 endif
310    
311     tp6 = MixingMap(i,j)%sigma6/rcut6
312     tp12 = tp6**2
313 gezelter 590 MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
314     MixingMap(i,j)%shiftedPot = doShift
315 gezelter 507
316 gezelter 590 if (i.ne.j) then
317     MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
318     MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
319     MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
320     MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
321     MixingMap(j,i)%delta = MixingMap(i,j)%delta
322     MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
323     MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
324     MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
325     endif
326    
327 gezelter 572 enddo
328     enddo
329    
330 chrisfen 217 haveMixingMap = .true.
331 gezelter 572
332 gezelter 135 end subroutine createMixingMap
333 gezelter 572
334 gezelter 115 subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
335     pot, f, do_pot)
336    
337     integer, intent(in) :: atom1, atom2
338 gezelter 572 integer :: atid1, atid2, ljt1, ljt2
339 gezelter 115 real( kind = dp ), intent(in) :: rij, r2
340     real( kind = dp ) :: pot, sw, vpair
341     real( kind = dp ), dimension(3,nLocal) :: f
342     real( kind = dp ), intent(in), dimension(3) :: d
343     real( kind = dp ), intent(inout), dimension(3) :: fpair
344     logical, intent(in) :: do_pot
345    
346     ! local Variables
347     real( kind = dp ) :: drdx, drdy, drdz
348     real( kind = dp ) :: fx, fy, fz
349     real( kind = dp ) :: pot_temp, dudr
350     real( kind = dp ) :: sigma6
351     real( kind = dp ) :: epsilon
352     real( kind = dp ) :: r6
353     real( kind = dp ) :: t6
354     real( kind = dp ) :: t12
355     real( kind = dp ) :: delta
356 gezelter 572 logical :: isSoftCore, shiftedPot
357 gezelter 135 integer :: id1, id2, localError
358 gezelter 115
359 gezelter 135 if (.not.haveMixingMap) then
360 gezelter 572 call createMixingMap()
361 gezelter 135 endif
362    
363 gezelter 115 ! Look up the correct parameters in the mixing matrix
364     #ifdef IS_MPI
365 gezelter 572 atid1 = atid_Row(atom1)
366     atid2 = atid_Col(atom2)
367 gezelter 115 #else
368 gezelter 572 atid1 = atid(atom1)
369     atid2 = atid(atom2)
370 gezelter 115 #endif
371    
372 gezelter 572 ljt1 = LJMap%atidToLJtype(atid1)
373     ljt2 = LJMap%atidToLJtype(atid2)
374    
375     sigma6 = MixingMap(ljt1,ljt2)%sigma6
376     epsilon = MixingMap(ljt1,ljt2)%epsilon
377     delta = MixingMap(ljt1,ljt2)%delta
378     isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
379     shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
380    
381 gezelter 115 r6 = r2 * r2 * r2
382 gezelter 507
383 gezelter 115 t6 = sigma6/ r6
384     t12 = t6 * t6
385 tim 378
386 gezelter 572 if (isSoftCore) then
387    
388     pot_temp = 4.0E0_DP * epsilon * t6
389     if (shiftedPot) then
390     pot_temp = pot_temp + delta
391     endif
392    
393     vpair = vpair + pot_temp
394    
395     dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
396    
397 tim 378 else
398 gezelter 507 pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
399 gezelter 572 if (shiftedPot) then
400 gezelter 507 pot_temp = pot_temp + delta
401     endif
402 gezelter 572
403 gezelter 507 vpair = vpair + pot_temp
404 gezelter 572
405 gezelter 507 dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
406 gezelter 115 endif
407 gezelter 507
408 gezelter 115 drdx = d(1) / rij
409     drdy = d(2) / rij
410     drdz = d(3) / rij
411 gezelter 507
412 gezelter 115 fx = dudr * drdx
413     fy = dudr * drdy
414     fz = dudr * drdz
415 gezelter 507
416 gezelter 115 #ifdef IS_MPI
417     if (do_pot) then
418 gezelter 662 pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
419     pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
420 gezelter 115 endif
421 gezelter 507
422 gezelter 115 f_Row(1,atom1) = f_Row(1,atom1) + fx
423     f_Row(2,atom1) = f_Row(2,atom1) + fy
424     f_Row(3,atom1) = f_Row(3,atom1) + fz
425 gezelter 507
426 gezelter 115 f_Col(1,atom2) = f_Col(1,atom2) - fx
427     f_Col(2,atom2) = f_Col(2,atom2) - fy
428     f_Col(3,atom2) = f_Col(3,atom2) - fz
429 gezelter 507
430 gezelter 115 #else
431     if (do_pot) pot = pot + sw*pot_temp
432    
433     f(1,atom1) = f(1,atom1) + fx
434     f(2,atom1) = f(2,atom1) + fy
435     f(3,atom1) = f(3,atom1) + fz
436 gezelter 507
437 gezelter 115 f(1,atom2) = f(1,atom2) - fx
438     f(2,atom2) = f(2,atom2) - fy
439     f(3,atom2) = f(3,atom2) - fz
440     #endif
441 gezelter 507
442 gezelter 115 #ifdef IS_MPI
443     id1 = AtomRowToGlobal(atom1)
444     id2 = AtomColToGlobal(atom2)
445     #else
446     id1 = atom1
447     id2 = atom2
448     #endif
449    
450     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
451 gezelter 507
452 gezelter 115 fpair(1) = fpair(1) + fx
453     fpair(2) = fpair(2) + fy
454     fpair(3) = fpair(3) + fz
455    
456     endif
457    
458     return
459 gezelter 507
460 gezelter 115 end subroutine do_lj_pair
461 gezelter 507
462 chuckv 492 subroutine destroyLJTypes()
463 gezelter 572
464     LJMap%nLJtypes = 0
465     LJMap%currentLJtype = 0
466    
467     if (associated(LJMap%LJtypes)) then
468     deallocate(LJMap%LJtypes)
469     LJMap%LJtypes => null()
470     end if
471    
472     if (associated(LJMap%atidToLJtype)) then
473     deallocate(LJMap%atidToLJtype)
474     LJMap%atidToLJtype => null()
475     end if
476    
477 chuckv 492 haveMixingMap = .false.
478     end subroutine destroyLJTypes
479    
480 gezelter 115 end module lj