ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90
Revision: 1930
Committed: Wed Jan 12 22:41:40 2005 UTC (19 years, 7 months ago) by gezelter
File size: 12107 byte(s)
Log Message:
merging new_design branch into OOPSE-2.0

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 1930 !! @version $Id: LJ.F90,v 1.6 2005-01-12 22:40:44 gezelter Exp $, $Date: 2005-01-12 22:40:44 $, $Name: not supported by cvs2svn $, $Revision: 1.6 $
47 gezelter 1608
48 gezelter 1930
49 gezelter 1608 module lj
50     use atype_module
51     use switcheroo
52     use vector_class
53     use simulation
54 gezelter 1628 use status
55 gezelter 1608 #ifdef IS_MPI
56     use mpiSimulation
57     #endif
58     use force_globals
59    
60     implicit none
61     PRIVATE
62 chuckv 1624
63     integer, parameter :: DP = selected_real_kind(15)
64 gezelter 1608
65 gezelter 1628 type, private :: LjType
66 gezelter 1930 integer :: c_ident
67     integer :: atid
68 gezelter 1628 real(kind=dp) :: sigma
69     real(kind=dp) :: epsilon
70     end type LjType
71 gezelter 1608
72 gezelter 1628 type(LjType), dimension(:), allocatable :: ParameterMap
73 gezelter 1608
74 gezelter 1628 logical, save :: haveMixingMap = .false.
75 gezelter 1608
76 gezelter 1628 type :: MixParameters
77     real(kind=DP) :: sigma
78     real(kind=DP) :: epsilon
79     real(kind=dp) :: sigma6
80     real(kind=dp) :: tp6
81     real(kind=dp) :: tp12
82     real(kind=dp) :: delta
83     end type MixParameters
84 chuckv 1624
85 gezelter 1628 type(MixParameters), dimension(:,:), allocatable :: MixingMap
86 chuckv 1624
87 gezelter 1628 real(kind=DP), save :: LJ_rcut
88     logical, save :: have_rcut = .false.
89     logical, save :: LJ_do_shift = .false.
90     logical, save :: useGeometricDistanceMixing = .false.
91    
92     !! Public methods and data
93 chuckv 1624
94 gezelter 1628 public :: setCutoffLJ
95     public :: useGeometricMixing
96     public :: do_lj_pair
97     public :: newLJtype
98 gezelter 1633 public :: getSigma
99     public :: getEpsilon
100 chuckv 1624
101 gezelter 1608 contains
102    
103 gezelter 1930 subroutine newLJtype(c_ident, sigma, epsilon, status)
104     integer,intent(in) :: c_ident
105 gezelter 1628 real(kind=dp),intent(in) :: sigma
106     real(kind=dp),intent(in) :: epsilon
107 chuckv 1624 integer,intent(out) :: status
108 gezelter 1930 integer :: nATypes, myATID
109     integer, pointer :: MatchList(:) => null()
110 gezelter 1628
111 chuckv 1624 status = 0
112    
113 gezelter 1930 myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
114 chuckv 1624
115 gezelter 1628 if (.not.allocated(ParameterMap)) then
116    
117 gezelter 1930 !call getMatchingElementList(atypes, "is_LennardJones", .true., &
118     ! nLJTypes, MatchList)
119 gezelter 1628 nAtypes = getSize(atypes)
120 chuckv 1624 if (nAtypes == 0) then
121 gezelter 1628 status = -1
122     return
123 chuckv 1624 end if
124 gezelter 1628
125     if (.not. allocated(ParameterMap)) then
126     allocate(ParameterMap(nAtypes))
127     endif
128    
129 chuckv 1624 end if
130    
131 gezelter 1930 if (myATID .gt. size(ParameterMap)) then
132 gezelter 1628 status = -1
133     return
134     endif
135 chuckv 1624
136 gezelter 1628 ! set the values for ParameterMap for this atom type:
137    
138 gezelter 1930 ParameterMap(myATID)%c_ident = c_ident
139     ParameterMap(myATID)%atid = myATID
140     ParameterMap(myATID)%epsilon = epsilon
141     ParameterMap(myATID)%sigma = sigma
142 chuckv 1624
143     end subroutine newLJtype
144 gezelter 1633
145     function getSigma(atid) result (s)
146     integer, intent(in) :: atid
147     integer :: localError
148     real(kind=dp) :: s
149 gezelter 1608
150 gezelter 1633 if (.not.allocated(ParameterMap)) then
151     call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
152     return
153     end if
154    
155     s = ParameterMap(atid)%sigma
156     end function getSigma
157    
158     function getEpsilon(atid) result (e)
159     integer, intent(in) :: atid
160     integer :: localError
161     real(kind=dp) :: e
162    
163     if (.not.allocated(ParameterMap)) then
164 gezelter 1930 call handleError("LJ", "no ParameterMap was present before first call of getEpsilon!")
165 gezelter 1633 return
166     end if
167    
168     e = ParameterMap(atid)%epsilon
169     end function getEpsilon
170    
171    
172 gezelter 1608 subroutine setCutoffLJ(rcut, do_shift, status)
173     logical, intent(in):: do_shift
174     integer :: status, myStatus
175     real(kind=dp) :: rcut
176    
177     #define __FORTRAN90
178     #include "UseTheForce/fSwitchingFunction.h"
179    
180     status = 0
181    
182     LJ_rcut = rcut
183     LJ_do_shift = do_shift
184     call set_switch(LJ_SWITCH, rcut, rcut)
185 gezelter 1628 have_rcut = .true.
186 gezelter 1608
187     return
188     end subroutine setCutoffLJ
189 gezelter 1628
190     subroutine useGeometricMixing()
191     useGeometricDistanceMixing = .true.
192     haveMixingMap = .false.
193     return
194     end subroutine useGeometricMixing
195 gezelter 1608
196 gezelter 1628 subroutine createMixingMap(status)
197 gezelter 1930 integer :: nATIDs
198 gezelter 1608 integer :: status
199     integer :: i
200     integer :: j
201 gezelter 1628 real ( kind = dp ) :: Sigma_i, Sigma_j
202     real ( kind = dp ) :: Epsilon_i, Epsilon_j
203 gezelter 1608 real ( kind = dp ) :: rcut6
204 gezelter 1628
205 gezelter 1608 status = 0
206    
207 gezelter 1930 nATIDs = size(ParameterMap)
208 gezelter 1628
209 gezelter 1930 if (nATIDs == 0) then
210 gezelter 1608 status = -1
211     return
212     end if
213 gezelter 1628
214     if (.not.have_rcut) then
215     status = -1
216     return
217 gezelter 1608 endif
218 gezelter 1628
219     if (.not. allocated(MixingMap)) then
220 gezelter 1930 allocate(MixingMap(nATIDs, nATIDs))
221 gezelter 1628 endif
222    
223 gezelter 1608 rcut6 = LJ_rcut**6
224 gezelter 1628
225     ! This loops through all atypes, even those that don't support LJ forces.
226 gezelter 1930 do i = 1, nATIDs
227 gezelter 1628
228     Epsilon_i = ParameterMap(i)%epsilon
229     Sigma_i = ParameterMap(i)%sigma
230    
231     ! do self mixing rule
232     MixingMap(i,i)%sigma = Sigma_i
233     MixingMap(i,i)%sigma6 = Sigma_i ** 6
234     MixingMap(i,i)%tp6 = (MixingMap(i,i)%sigma6)/rcut6
235     MixingMap(i,i)%tp12 = (MixingMap(i,i)%tp6) ** 2
236     MixingMap(i,i)%epsilon = Epsilon_i
237     MixingMap(i,i)%delta = -4.0_DP * MixingMap(i,i)%epsilon * &
238     (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
239    
240 gezelter 1930 do j = i + 1, nATIDs
241 chuckv 1624
242 gezelter 1628 Epsilon_j = ParameterMap(j)%epsilon
243     Sigma_j = ParameterMap(j)%sigma
244 gezelter 1608
245 gezelter 1628 ! only the distance parameter uses different mixing policies
246     if (useGeometricDistanceMixing) then
247     ! only for OPLS as far as we can tell
248     MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
249     else
250     ! everyone else
251     MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
252     endif
253 gezelter 1608
254 gezelter 1628 ! energy parameter is always geometric mean:
255     MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
256    
257     MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
258     MixingMap(i,j)%tp6 = MixingMap(i,j)%sigma6/rcut6
259     MixingMap(i,j)%tp12 = (MixingMap(i,j)%tp6) ** 2
260 gezelter 1608
261 gezelter 1628 MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
262     (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
263 gezelter 1608
264 gezelter 1628 MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
265     MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
266     MixingMap(j,i)%tp6 = MixingMap(i,j)%tp6
267     MixingMap(j,i)%tp12 = MixingMap(i,j)%tp12
268     MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
269     MixingMap(j,i)%delta = MixingMap(i,j)%delta
270 gezelter 1608
271 gezelter 1628 end do
272 gezelter 1608 end do
273    
274 chrisfen 1754 haveMixingMap = .true.
275    
276 gezelter 1628 end subroutine createMixingMap
277    
278 gezelter 1608 subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
279     pot, f, do_pot)
280    
281     integer, intent(in) :: atom1, atom2
282     real( kind = dp ), intent(in) :: rij, r2
283     real( kind = dp ) :: pot, sw, vpair
284     real( kind = dp ), dimension(3,nLocal) :: f
285     real( kind = dp ), intent(in), dimension(3) :: d
286     real( kind = dp ), intent(inout), dimension(3) :: fpair
287     logical, intent(in) :: do_pot
288    
289     ! local Variables
290     real( kind = dp ) :: drdx, drdy, drdz
291     real( kind = dp ) :: fx, fy, fz
292     real( kind = dp ) :: pot_temp, dudr
293     real( kind = dp ) :: sigma6
294     real( kind = dp ) :: epsilon
295     real( kind = dp ) :: r6
296     real( kind = dp ) :: t6
297     real( kind = dp ) :: t12
298     real( kind = dp ) :: delta
299 gezelter 1628 integer :: id1, id2, localError
300 gezelter 1608
301 gezelter 1628 if (.not.haveMixingMap) then
302     localError = 0
303     call createMixingMap(localError)
304     if ( localError .ne. 0 ) then
305     call handleError("LJ", "MixingMap creation failed!")
306     return
307     end if
308     endif
309    
310 gezelter 1608 ! Look up the correct parameters in the mixing matrix
311     #ifdef IS_MPI
312 gezelter 1628 sigma6 = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
313     epsilon = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
314     delta = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
315 gezelter 1608 #else
316 gezelter 1628 sigma6 = MixingMap(atid(atom1),atid(atom2))%sigma6
317     epsilon = MixingMap(atid(atom1),atid(atom2))%epsilon
318     delta = MixingMap(atid(atom1),atid(atom2))%delta
319 gezelter 1608 #endif
320    
321     r6 = r2 * r2 * r2
322    
323     t6 = sigma6/ r6
324     t12 = t6 * t6
325    
326     pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
327     if (LJ_do_shift) then
328     pot_temp = pot_temp + delta
329     endif
330    
331     vpair = vpair + pot_temp
332    
333     dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
334    
335     drdx = d(1) / rij
336     drdy = d(2) / rij
337     drdz = d(3) / rij
338    
339     fx = dudr * drdx
340     fy = dudr * drdy
341     fz = dudr * drdz
342    
343    
344     #ifdef IS_MPI
345     if (do_pot) then
346     pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
347     pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
348     endif
349    
350     f_Row(1,atom1) = f_Row(1,atom1) + fx
351     f_Row(2,atom1) = f_Row(2,atom1) + fy
352     f_Row(3,atom1) = f_Row(3,atom1) + fz
353    
354     f_Col(1,atom2) = f_Col(1,atom2) - fx
355     f_Col(2,atom2) = f_Col(2,atom2) - fy
356     f_Col(3,atom2) = f_Col(3,atom2) - fz
357    
358     #else
359     if (do_pot) pot = pot + sw*pot_temp
360    
361     f(1,atom1) = f(1,atom1) + fx
362     f(2,atom1) = f(2,atom1) + fy
363     f(3,atom1) = f(3,atom1) + fz
364    
365     f(1,atom2) = f(1,atom2) - fx
366     f(2,atom2) = f(2,atom2) - fy
367     f(3,atom2) = f(3,atom2) - fz
368     #endif
369    
370     #ifdef IS_MPI
371     id1 = AtomRowToGlobal(atom1)
372     id2 = AtomColToGlobal(atom2)
373     #else
374     id1 = atom1
375     id2 = atom2
376     #endif
377    
378     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
379    
380     fpair(1) = fpair(1) + fx
381     fpair(2) = fpair(2) + fy
382     fpair(3) = fpair(3) + fz
383    
384     endif
385    
386     return
387    
388     end subroutine do_lj_pair
389    
390    
391     !! Calculates the mixing for sigma or epslon
392    
393     end module lj
394 chuckv 1624
395 gezelter 1930 subroutine newLJtype(c_ident, sigma, epsilon, status)
396 gezelter 1628 use lj, ONLY : module_newLJtype => newLJtype
397     integer, parameter :: DP = selected_real_kind(15)
398 gezelter 1930 integer,intent(inout) :: c_ident
399 gezelter 1628 real(kind=dp),intent(inout) :: sigma
400     real(kind=dp),intent(inout) :: epsilon
401     integer,intent(inout) :: status
402    
403 gezelter 1930 call module_newLJtype(c_ident, sigma, epsilon, status)
404 gezelter 1628
405     end subroutine newLJtype
406 chuckv 1624
407 gezelter 1628 subroutine useGeometricMixing()
408     use lj, ONLY: module_useGeometricMixing => useGeometricMixing
409    
410     call module_useGeometricMixing()
411     return
412     end subroutine useGeometricMixing