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 2062 by tim, Fri Feb 25 21:22:00 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.8 2005-02-25 21:22:00 tim Exp $, $Date: 2005-02-25 21:22:00 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
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 >     logical :: soft_pot
71 >  end type LjType
72    
73 +  type(LjType), dimension(:), allocatable :: ParameterMap
74 +  
75 +  logical, save :: haveMixingMap = .false.
76 +  
77 +  type :: MixParameters
78 +     real(kind=DP) :: sigma
79 +     real(kind=DP) :: epsilon
80 +     real(kind=dp)  :: sigma6
81 +     real(kind=dp)  :: tp6
82 +     real(kind=dp)  :: tp12
83 +     real(kind=dp)  :: delta
84 +     logical :: soft_pot    
85 +  end type MixParameters
86 +  
87 +  type(MixParameters), dimension(:,:), allocatable :: MixingMap
88 +  
89 +  real(kind=DP), save :: LJ_rcut
90 +  logical, save :: have_rcut = .false.
91 +  logical, save :: LJ_do_shift = .false.
92 +  logical, save :: useGeometricDistanceMixing = .false.
93 +  
94    !! Public methods and data
95 <  public :: init_LJ_FF
95 >  
96    public :: setCutoffLJ
97 +  public :: useGeometricMixing
98    public :: do_lj_pair
99 +  public :: newLJtype  
100 +  public :: getSigma
101 +  public :: getEpsilon
102    
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 –  
103   contains
104  
105 <  subroutine init_LJ_FF(mix_Policy, status)
106 <    integer, intent(in) :: mix_Policy
107 <    integer, intent(out) :: status
108 <    integer :: myStatus
105 >  subroutine newLJtype(c_ident, sigma, epsilon, soft_pot, status)
106 >    integer,intent(in) :: c_ident
107 >    real(kind=dp),intent(in) :: sigma
108 >    real(kind=dp),intent(in) :: epsilon
109 >    integer, intent(in) :: soft_pot
110 >    integer,intent(out) :: status
111 >    integer :: nATypes, myATID
112 >    integer, pointer :: MatchList(:) => null()
113 >
114 >    status = 0
115      
116 <    if (mix_Policy == LB_MIXING_RULE) then
117 <       LJ_Mixing_Policy = LB_MIXING_RULE
118 <    else
119 <       if (mix_Policy == EXPLICIT_MIXING_RULE) then
120 <          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
121 <       else
122 <          write(*,*) 'Unknown Mixing Policy!'
116 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
117 >
118 >    if (.not.allocated(ParameterMap)) then
119 >      
120 >       !call getMatchingElementList(atypes, "is_LennardJones", .true., &
121 >       !     nLJTypes, MatchList)
122 >       nAtypes = getSize(atypes)
123 >       if (nAtypes == 0) then
124            status = -1
125            return
126 +       end if
127 +      
128 +       if (.not. allocated(ParameterMap)) then
129 +          allocate(ParameterMap(nAtypes))
130         endif
131 +      
132 +    end if
133 +
134 +    if (myATID .gt. size(ParameterMap)) then
135 +       status = -1
136 +       return
137      endif
138 +    
139 +    ! set the values for ParameterMap for this atom type:
140  
141 <    havePolicy = .true.
141 >    ParameterMap(myATID)%c_ident = c_ident
142 >    ParameterMap(myATID)%atid = myATID
143 >    ParameterMap(myATID)%epsilon = epsilon
144 >    ParameterMap(myATID)%sigma = sigma
145 >    if (soft_pot == 1) then
146 >      ParameterMap(myATID)%soft_pot = .true.
147 >    else
148 >      ParameterMap(myATID)%soft_pot = .false.
149 >    endif
150 >  end subroutine newLJtype
151  
152 <    if (haveCut) then
153 <       status = 0
154 <       call createMixingList(myStatus)
155 <       if (myStatus /= 0) then
156 <          status = -1
157 <          return
158 <       end if
159 <      
83 <       LJ_FF_initialized = .true.
152 >  function getSigma(atid) result (s)
153 >    integer, intent(in) :: atid
154 >    integer :: localError
155 >    real(kind=dp) :: s
156 >    
157 >    if (.not.allocated(ParameterMap)) then
158 >       call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
159 >       return
160      end if
161 <  
162 <  end subroutine init_LJ_FF
163 <  
161 >    
162 >    s = ParameterMap(atid)%sigma
163 >  end function getSigma
164 >
165 >  function getEpsilon(atid) result (e)
166 >    integer, intent(in) :: atid
167 >    integer :: localError
168 >    real(kind=dp) :: e
169 >    
170 >    if (.not.allocated(ParameterMap)) then
171 >       call handleError("LJ", "no ParameterMap was present before first call of getEpsilon!")
172 >       return
173 >    end if
174 >    
175 >    e = ParameterMap(atid)%epsilon
176 >  end function getEpsilon
177 >
178 >
179    subroutine setCutoffLJ(rcut, do_shift, status)
180      logical, intent(in):: do_shift
181      integer :: status, myStatus
# Line 98 | Line 189 | contains
189      LJ_rcut = rcut
190      LJ_do_shift = do_shift
191      call set_switch(LJ_SWITCH, rcut, rcut)
192 <    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    
192 >    have_rcut = .true.
193      
194      return
195    end subroutine setCutoffLJ
196 +
197 +  subroutine useGeometricMixing()
198 +    useGeometricDistanceMixing = .true.
199 +    haveMixingMap = .false.
200 +    return
201 +  end subroutine useGeometricMixing
202    
203 <  subroutine createMixingList(status)
204 <    integer :: nAtypes
203 >  subroutine createMixingMap(status)
204 >    integer :: nATIDs
205      integer :: status
206      integer :: i
207      integer :: j
208 <    real ( kind = dp ) :: mySigma_i,mySigma_j
209 <    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
208 >    real ( kind = dp ) :: Sigma_i, Sigma_j
209 >    real ( kind = dp ) :: Epsilon_i, Epsilon_j
210      real ( kind = dp ) :: rcut6
211 <    logical :: I_isLJ, J_isLJ
211 >
212      status = 0
213      
214 <    nAtypes = getSize(atypes)
215 <    if (nAtypes == 0) then
214 >    nATIDs = size(ParameterMap)
215 >    
216 >    if (nATIDs == 0) then
217         status = -1
218         return
219      end if
133 –        
134 –    if (.not. associated(ljMixed)) then
135 –       allocate(ljMixed(nAtypes, nAtypes))
136 –    endif
220  
221 +    if (.not.have_rcut) then
222 +       status = -1
223 +       return
224 +    endif
225 +    
226 +    if (.not. allocated(MixingMap)) then
227 +       allocate(MixingMap(nATIDs, nATIDs))
228 +    endif
229 +    
230      rcut6 = LJ_rcut**6
231 <
232 < ! This loops through all atypes, even those that don't support LJ forces.
233 <    do i = 1, nAtypes
234 <
235 <       call getElementProperty(atypes, i, "is_LJ", I_isLJ)
236 <
237 <       if (I_isLJ) then
238 <
239 <          call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
240 <          call getElementProperty(atypes, i, "lj_sigma",   mySigma_i)
241 <          ! do self mixing rule
242 <          ljMixed(i,i)%sigma   = mySigma_i
231 >    
232 >    ! This loops through all atypes, even those that don't support LJ forces.
233 >    do i = 1, nATIDs
234 >      
235 >       Epsilon_i = ParameterMap(i)%epsilon
236 >       Sigma_i = ParameterMap(i)%sigma
237 >      
238 >       ! do self mixing rule
239 >       MixingMap(i,i)%sigma   = Sigma_i          
240 >       MixingMap(i,i)%sigma6  = Sigma_i ** 6          
241 >       MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
242 >       MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
243 >       MixingMap(i,i)%epsilon = Epsilon_i          
244 >       MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
245 >            (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
246 >       MixingMap(i,i)%soft_pot = ParameterMap(i)%soft_pot
247 >      
248 >       do j = i + 1, nATIDs
249            
250 <          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
250 >          Epsilon_j = ParameterMap(j)%epsilon
251 >          Sigma_j = ParameterMap(j)%sigma
252            
253 <          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
253 >          ! only the distance parameter uses different mixing policies
254 >          if (useGeometricDistanceMixing) then
255 >             ! only for OPLS as far as we can tell
256 >             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
257 >          else
258 >             ! everyone else
259 >             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
260 >          endif
261            
262 <          ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
262 >          ! energy parameter is always geometric mean:
263 >          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
264 >                    
265 >          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
266 >          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
267 >          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
268            
269 +          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
270 +               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
271 +
272 +          MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
273 +
274            
275 <          ljMixed(i,i)%epsilon = myEpsilon_i
275 >          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
276 >          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
277 >          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
278 >          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
279 >          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
280 >          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
281 >          MixingMap(j,i)%soft_pot   = MixingMap(i,j)%soft_pot
282            
283 <          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
283 >       end do
284      end do
285      
286 <  end subroutine createMixingList
287 <  
286 >    haveMixingMap = .true.
287 >
288 >  end subroutine createMixingMap
289 >        
290    subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
291         pot, f, do_pot)
292  
# Line 225 | Line 308 | contains
308      real( kind = dp ) :: t6
309      real( kind = dp ) :: t12
310      real( kind = dp ) :: delta
311 <    integer :: id1, id2
311 >    logical :: soft_pot
312 >    integer :: id1, id2, localError
313  
314 +    if (.not.haveMixingMap) then
315 +       localError = 0
316 +       call createMixingMap(localError)
317 +       if ( localError .ne. 0 ) then
318 +          call handleError("LJ", "MixingMap creation failed!")
319 +          return
320 +       end if
321 +    endif
322 +
323      ! Look up the correct parameters in the mixing matrix
324   #ifdef IS_MPI
325 <    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
326 <    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
327 <    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
325 >    sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
326 >    epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
327 >    delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
328 >    soft_pot =  MixingMap(atid_Row(atom1),atid_Col(atom2))%soft_pot
329   #else
330 <    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
331 <    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
332 <    delta    = ljMixed(atid(atom1),atid(atom2))%delta
330 >    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
331 >    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
332 >    delta    = MixingMap(atid(atom1),atid(atom2))%delta
333 >    soft_pot    = MixingMap(atid(atom1),atid(atom2))%soft_pot
334   #endif
335  
336      r6 = r2 * r2 * r2
337      
338      t6  = sigma6/ r6
339      t12 = t6 * t6    
245 –  
246 –    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
247 –    if (LJ_do_shift) then
248 –       pot_temp = pot_temp + delta
249 –    endif
340  
341 <    vpair = vpair + pot_temp
341 >    if (soft_pot) then
342 >      pot_temp = 4.0E0_DP * epsilon * t12
343 >      if (LJ_do_shift) then
344 >         pot_temp = pot_temp + delta
345 >      endif
346 >  
347 >      vpair = vpair + pot_temp
348 >        
349 >      dudr = sw * 24.0E0_DP * epsilon * (-2.0E0_DP)*t12 / rij
350 >
351 >    else
352 >      pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
353 >      if (LJ_do_shift) then
354 >         pot_temp = pot_temp + delta
355 >      endif
356 >  
357 >      vpair = vpair + pot_temp
358 >        
359 >      dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
360 >    endif
361        
253 –    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
254 –      
362      drdx = d(1) / rij
363      drdy = d(2) / rij
364      drdz = d(3) / rij
# Line 309 | Line 416 | contains
416    
417    
418    !! Calculates the mixing for sigma or epslon
312 –  
313 –  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
314 –    character(len=*) :: thisParam
315 –    real(kind = dp)  :: param1
316 –    real(kind = dp)  :: param2
317 –    real(kind = dp ) :: myMixParam
318 –
319 –    integer, optional :: status  
320 –
321 –    myMixParam = 0.0_dp
419      
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
338 –  
420   end module lj

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines