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

Diff Legend

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