ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90 (file contents):
Revision 1608 by gezelter, Wed Oct 20 04:02:48 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC

# Line 1 | Line 1
1 + !!
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   !! Calculates Long Range forces Lennard-Jones interactions.
2 !! Corresponds to the force field defined in lj_FF.cpp
44   !! @author Charles F. Vardeman II
45   !! @author Matthew Meineke
46 < !! @version $Id: LJ.F90,v 1.1 2004-10-20 04:02:48 gezelter Exp $, $Date: 2004-10-20 04:02:48 $, $Name: not supported by cvs2svn $, $Revision: 1.1 $
46 > !! @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  
48 +
49   module lj
8  use definitions
50    use atype_module
51    use switcheroo
52    use vector_class
53    use simulation
54 +  use status
55   #ifdef IS_MPI
56    use mpiSimulation
57   #endif
# Line 17 | Line 59 | module lj
59  
60    implicit none
61    PRIVATE
20
21 #define __FORTRAN90
22 #include "UseTheForce/fForceField.h"
23
24  integer, save :: LJ_Mixing_Policy
25  real(kind=DP), save :: LJ_rcut
26  logical, save :: havePolicy = .false.
27  logical, save :: haveCut = .false.
28  logical, save :: LJ_do_shift = .false.
62    
63 <  !! Logical has lj force field module been initialized?
63 >  integer, parameter :: DP = selected_real_kind(15)
64    
65 <  logical, save :: LJ_FF_initialized = .false.
65 >  type, private :: LjType
66 >     integer :: c_ident
67 >     integer :: atid
68 >     real(kind=dp) :: sigma
69 >     real(kind=dp) :: epsilon
70 >  end type LjType
71    
72 +  type(LjType), dimension(:), allocatable :: ParameterMap
73 +  
74 +  logical, save :: haveMixingMap = .false.
75 +  
76 +  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 +  
85 +  type(MixParameters), dimension(:,:), allocatable :: MixingMap
86 +  
87 +  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 <  public :: init_LJ_FF
93 >  
94    public :: setCutoffLJ
95 +  public :: useGeometricMixing
96    public :: do_lj_pair
97 +  public :: newLJtype  
98 +  public :: getSigma
99 +  public :: getEpsilon
100    
39  type :: lj_mixed_params
40     !! Lennard-Jones epsilon
41     real ( kind = dp )  :: epsilon = 0.0_dp
42     !! Lennard-Jones Sigma
43     real ( kind = dp )  :: sigma = 0.0_dp
44     !! Lennard-Jones Sigma to sixth
45     real ( kind = dp )  :: sigma6 = 0.0_dp
46     !!
47     real ( kind = dp )  :: tp6
48     real ( kind = dp )  :: tp12
49     real ( kind = dp )  :: delta  = 0.0_dp
50  end type lj_mixed_params
51  
52  type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
53  
101   contains
102  
103 <  subroutine init_LJ_FF(mix_Policy, status)
104 <    integer, intent(in) :: mix_Policy
105 <    integer, intent(out) :: status
106 <    integer :: myStatus
103 >  subroutine newLJtype(c_ident, sigma, epsilon, status)
104 >    integer,intent(in) :: c_ident
105 >    real(kind=dp),intent(in) :: sigma
106 >    real(kind=dp),intent(in) :: epsilon
107 >    integer,intent(out) :: status
108 >    integer :: nATypes, myATID
109 >    integer, pointer :: MatchList(:) => null()
110 >
111 >    status = 0
112      
113 <    if (mix_Policy == LB_MIXING_RULE) then
114 <       LJ_Mixing_Policy = LB_MIXING_RULE
115 <    else
116 <       if (mix_Policy == EXPLICIT_MIXING_RULE) then
117 <          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
118 <       else
119 <          write(*,*) 'Unknown Mixing Policy!'
113 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
114 >
115 >    if (.not.allocated(ParameterMap)) then
116 >      
117 >       !call getMatchingElementList(atypes, "is_LennardJones", .true., &
118 >       !     nLJTypes, MatchList)
119 >       nAtypes = getSize(atypes)
120 >       if (nAtypes == 0) then
121            status = -1
122            return
123 +       end if
124 +      
125 +       if (.not. allocated(ParameterMap)) then
126 +          allocate(ParameterMap(nAtypes))
127         endif
128 +      
129 +    end if
130 +
131 +    if (myATID .gt. size(ParameterMap)) then
132 +       status = -1
133 +       return
134      endif
135 +    
136 +    ! set the values for ParameterMap for this atom type:
137  
138 <    havePolicy = .true.
138 >    ParameterMap(myATID)%c_ident = c_ident
139 >    ParameterMap(myATID)%atid = myATID
140 >    ParameterMap(myATID)%epsilon = epsilon
141 >    ParameterMap(myATID)%sigma = sigma
142 >    
143 >  end subroutine newLJtype
144  
145 <    if (haveCut) then
146 <       status = 0
147 <       call createMixingList(myStatus)
148 <       if (myStatus /= 0) then
149 <          status = -1
150 <          return
151 <       end if
152 <      
83 <       LJ_FF_initialized = .true.
145 >  function getSigma(atid) result (s)
146 >    integer, intent(in) :: atid
147 >    integer :: localError
148 >    real(kind=dp) :: s
149 >    
150 >    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 <  end subroutine init_LJ_FF
156 <  
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 >       call handleError("LJ", "no ParameterMap was present before first call of getEpsilon!")
165 >       return
166 >    end if
167 >    
168 >    e = ParameterMap(atid)%epsilon
169 >  end function getEpsilon
170 >
171 >
172    subroutine setCutoffLJ(rcut, do_shift, status)
173      logical, intent(in):: do_shift
174      integer :: status, myStatus
# Line 98 | Line 182 | contains
182      LJ_rcut = rcut
183      LJ_do_shift = do_shift
184      call set_switch(LJ_SWITCH, rcut, rcut)
185 <    haveCut = .true.
102 <
103 <    if (havePolicy) then
104 <       status = 0
105 <       call createMixingList(myStatus)
106 <       if (myStatus /= 0) then
107 <          status = -1
108 <          return
109 <       end if
110 <      
111 <       LJ_FF_initialized = .true.
112 <    end if    
185 >    have_rcut = .true.
186      
187      return
188    end subroutine setCutoffLJ
189 +
190 +  subroutine useGeometricMixing()
191 +    useGeometricDistanceMixing = .true.
192 +    haveMixingMap = .false.
193 +    return
194 +  end subroutine useGeometricMixing
195    
196 <  subroutine createMixingList(status)
197 <    integer :: nAtypes
196 >  subroutine createMixingMap(status)
197 >    integer :: nATIDs
198      integer :: status
199      integer :: i
200      integer :: j
201 <    real ( kind = dp ) :: mySigma_i,mySigma_j
202 <    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
201 >    real ( kind = dp ) :: Sigma_i, Sigma_j
202 >    real ( kind = dp ) :: Epsilon_i, Epsilon_j
203      real ( kind = dp ) :: rcut6
204 <    logical :: I_isLJ, J_isLJ
204 >
205      status = 0
206      
207 <    nAtypes = getSize(atypes)
208 <    if (nAtypes == 0) then
207 >    nATIDs = size(ParameterMap)
208 >    
209 >    if (nATIDs == 0) then
210         status = -1
211         return
212      end if
133        
134    if (.not. associated(ljMixed)) then
135       allocate(ljMixed(nAtypes, nAtypes))
136    endif
213  
214 +    if (.not.have_rcut) then
215 +       status = -1
216 +       return
217 +    endif
218 +    
219 +    if (.not. allocated(MixingMap)) then
220 +       allocate(MixingMap(nATIDs, nATIDs))
221 +    endif
222 +    
223      rcut6 = LJ_rcut**6
224 <
225 < ! This loops through all atypes, even those that don't support LJ forces.
226 <    do i = 1, nAtypes
227 <
228 <       call getElementProperty(atypes, i, "is_LJ", I_isLJ)
229 <
230 <       if (I_isLJ) then
231 <
232 <          call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
233 <          call getElementProperty(atypes, i, "lj_sigma",   mySigma_i)
234 <          ! do self mixing rule
235 <          ljMixed(i,i)%sigma   = mySigma_i
224 >    
225 >    ! This loops through all atypes, even those that don't support LJ forces.
226 >    do i = 1, nATIDs
227 >      
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 >       do j = i + 1, nATIDs
241            
242 <          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
242 >          Epsilon_j = ParameterMap(j)%epsilon
243 >          Sigma_j = ParameterMap(j)%sigma
244            
245 <          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
245 >          ! 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            
254 <          ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
254 >          ! 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            
261 +          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
262 +               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
263            
264 <          ljMixed(i,i)%epsilon = myEpsilon_i
264 >          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            
271 <          ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
162 <            (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
163 <          
164 <          do j = i + 1, nAtypes
165 <
166 <             call getElementProperty(atypes, j, "is_LJ", J_isLJ)
167 <            
168 <             if (J_isLJ) then                
169 <            
170 <                call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
171 <                call getElementProperty(atypes,j,"lj_sigma",  mySigma_j)
172 <                
173 <                ljMixed(i,j)%sigma  =  &
174 <                     calcLJMix("sigma",mySigma_i, &
175 <                     mySigma_j)
176 <                
177 <                ljMixed(i,j)%sigma6 = &
178 <                     (ljMixed(i,j)%sigma)**6
179 <                
180 <                
181 <                ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
182 <                
183 <                ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
184 <                
185 <                
186 <                ljMixed(i,j)%epsilon = &
187 <                     calcLJMix("epsilon",myEpsilon_i, &
188 <                     myEpsilon_j)
189 <                
190 <                ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
191 <                     (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
192 <                
193 <                
194 <                ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
195 <                ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
196 <                ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
197 <                ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
198 <                ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
199 <                ljMixed(j,i)%delta   = ljMixed(i,j)%delta
200 <             endif                
201 <          end do
202 <       endif
271 >       end do
272      end do
273      
274 <  end subroutine createMixingList
275 <  
274 >    haveMixingMap = .true.
275 >
276 >  end subroutine createMixingMap
277 >        
278    subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
279         pot, f, do_pot)
280  
# Line 225 | Line 296 | contains
296      real( kind = dp ) :: t6
297      real( kind = dp ) :: t12
298      real( kind = dp ) :: delta
299 <    integer :: id1, id2
299 >    integer :: id1, id2, localError
300  
301 +    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      ! Look up the correct parameters in the mixing matrix
311   #ifdef IS_MPI
312 <    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
313 <    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
314 <    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
312 >    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   #else
316 <    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
317 <    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
318 <    delta    = ljMixed(atid(atom1),atid(atom2))%delta
316 >    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
317 >    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
318 >    delta    = MixingMap(atid(atom1),atid(atom2))%delta
319   #endif
320  
321      r6 = r2 * r2 * r2
# Line 309 | Line 389 | contains
389    
390    
391    !! Calculates the mixing for sigma or epslon
392 +    
393 + end module lj
394 +
395 + subroutine newLJtype(c_ident, sigma, epsilon, status)
396 +  use lj, ONLY : module_newLJtype => newLJtype
397 +  integer, parameter :: DP = selected_real_kind(15)
398 +  integer,intent(inout) :: c_ident
399 +  real(kind=dp),intent(inout) :: sigma
400 +  real(kind=dp),intent(inout) :: epsilon
401 +  integer,intent(inout) :: status
402    
403 <  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
404 <    character(len=*) :: thisParam
405 <    real(kind = dp)  :: param1
316 <    real(kind = dp)  :: param2
317 <    real(kind = dp ) :: myMixParam
403 >  call module_newLJtype(c_ident, sigma, epsilon, status)
404 >  
405 > end subroutine newLJtype
406  
407 <    integer, optional :: status  
408 <
321 <    myMixParam = 0.0_dp
322 <    
323 <    if (present(status)) status = 0
324 <    select case (LJ_Mixing_Policy)
325 <    case (1)
326 <       select case (thisParam)
327 <       case ("sigma")
328 <          myMixParam = 0.5_dp * (param1 + param2)
329 <       case ("epsilon")
330 <          myMixParam = sqrt(param1 * param2)
331 <       case default
332 <          status = -1
333 <       end select
334 <    case default
335 <       status = -1
336 <    end select
337 <  end function calcLJMix
407 > subroutine useGeometricMixing()
408 >  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
409    
410 < end module lj
410 >  call module_useGeometricMixing()
411 >  return
412 > end subroutine useGeometricMixing

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines