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 1633 by gezelter, Fri Oct 22 20:22:48 2004 UTC

# Line 1 | Line 1
1   !! Calculates Long Range forces Lennard-Jones interactions.
2 !! Corresponds to the force field defined in lj_FF.cpp
2   !! @author Charles F. Vardeman II
3   !! @author Matthew Meineke
4 < !! @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 $
4 > !! @version $Id: LJ.F90,v 1.4 2004-10-22 20:22:47 gezelter Exp $, $Date: 2004-10-22 20:22:47 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
5  
6   module lj
8  use definitions
7    use atype_module
8    use switcheroo
9    use vector_class
10    use simulation
11 +  use status
12   #ifdef IS_MPI
13    use mpiSimulation
14   #endif
# Line 17 | Line 16 | module lj
16  
17    implicit none
18    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.
19    
20 <  !! Logical has lj force field module been initialized?
20 >  integer, parameter :: DP = selected_real_kind(15)
21    
22 <  logical, save :: LJ_FF_initialized = .false.
22 >  type, private :: LjType
23 >     integer :: ident
24 >     real(kind=dp) :: sigma
25 >     real(kind=dp) :: epsilon
26 >  end type LjType
27    
28 +  type(LjType), dimension(:), allocatable :: ParameterMap
29 +  
30 +  logical, save :: haveMixingMap = .false.
31 +  
32 +  type :: MixParameters
33 +     real(kind=DP) :: sigma
34 +     real(kind=DP) :: epsilon
35 +     real(kind=dp)  :: sigma6
36 +     real(kind=dp)  :: tp6
37 +     real(kind=dp)  :: tp12
38 +     real(kind=dp)  :: delta
39 +  end type MixParameters
40 +  
41 +  type(MixParameters), dimension(:,:), allocatable :: MixingMap
42 +  
43 +  real(kind=DP), save :: LJ_rcut
44 +  logical, save :: have_rcut = .false.
45 +  logical, save :: LJ_do_shift = .false.
46 +  logical, save :: useGeometricDistanceMixing = .false.
47 +  
48    !! Public methods and data
49 <  public :: init_LJ_FF
49 >  
50    public :: setCutoffLJ
51 +  public :: useGeometricMixing
52    public :: do_lj_pair
53 +  public :: newLJtype  
54 +  public :: getSigma
55 +  public :: getEpsilon
56    
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  
57   contains
58  
59 <  subroutine init_LJ_FF(mix_Policy, status)
60 <    integer, intent(in) :: mix_Policy
61 <    integer, intent(out) :: status
62 <    integer :: myStatus
59 >  subroutine newLJtype(ident, sigma, epsilon, status)
60 >    integer,intent(in) :: ident
61 >    real(kind=dp),intent(in) :: sigma
62 >    real(kind=dp),intent(in) :: epsilon
63 >    integer,intent(out) :: status
64 >    integer :: nAtypes
65 >
66 >    status = 0
67      
68 <    if (mix_Policy == LB_MIXING_RULE) then
69 <       LJ_Mixing_Policy = LB_MIXING_RULE
70 <    else
71 <       if (mix_Policy == EXPLICIT_MIXING_RULE) then
72 <          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
73 <       else
74 <          write(*,*) 'Unknown Mixing Policy!'
68 >    !! Be simple-minded and assume that we need a ParameterMap that
69 >    !! is the same size as the total number of atom types
70 >
71 >    if (.not.allocated(ParameterMap)) then
72 >      
73 >       nAtypes = getSize(atypes)
74 >    
75 >       if (nAtypes == 0) then
76            status = -1
77            return
78 +       end if
79 +      
80 +       if (.not. allocated(ParameterMap)) then
81 +          allocate(ParameterMap(nAtypes))
82         endif
83 +      
84 +    end if
85 +
86 +    if (ident .gt. size(ParameterMap)) then
87 +       status = -1
88 +       return
89      endif
90 +    
91 +    ! set the values for ParameterMap for this atom type:
92  
93 <    havePolicy = .true.
93 >    ParameterMap(ident)%ident = ident
94 >    ParameterMap(ident)%epsilon = epsilon
95 >    ParameterMap(ident)%sigma = sigma
96 >    
97 >  end subroutine newLJtype
98  
99 <    if (haveCut) then
100 <       status = 0
101 <       call createMixingList(myStatus)
102 <       if (myStatus /= 0) then
103 <          status = -1
104 <          return
105 <       end if
106 <      
83 <       LJ_FF_initialized = .true.
99 >  function getSigma(atid) result (s)
100 >    integer, intent(in) :: atid
101 >    integer :: localError
102 >    real(kind=dp) :: s
103 >    
104 >    if (.not.allocated(ParameterMap)) then
105 >       call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
106 >       return
107      end if
108 <  
109 <  end subroutine init_LJ_FF
110 <  
108 >    
109 >    s = ParameterMap(atid)%sigma
110 >  end function getSigma
111 >
112 >  function getEpsilon(atid) result (e)
113 >    integer, intent(in) :: atid
114 >    integer :: localError
115 >    real(kind=dp) :: e
116 >    
117 >    if (.not.allocated(ParameterMap)) then
118 >       call handleError("dipole-dipole", "no ParameterMap was present before first call of getEpsilon!")
119 >       return
120 >    end if
121 >    
122 >    e = ParameterMap(atid)%epsilon
123 >  end function getEpsilon
124 >
125 >
126    subroutine setCutoffLJ(rcut, do_shift, status)
127      logical, intent(in):: do_shift
128      integer :: status, myStatus
# Line 98 | Line 136 | contains
136      LJ_rcut = rcut
137      LJ_do_shift = do_shift
138      call set_switch(LJ_SWITCH, rcut, rcut)
139 <    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    
139 >    have_rcut = .true.
140      
141      return
142    end subroutine setCutoffLJ
143 +
144 +  subroutine useGeometricMixing()
145 +    useGeometricDistanceMixing = .true.
146 +    haveMixingMap = .false.
147 +    return
148 +  end subroutine useGeometricMixing
149    
150 <  subroutine createMixingList(status)
150 >  subroutine createMixingMap(status)
151      integer :: nAtypes
152      integer :: status
153      integer :: i
154      integer :: j
155 <    real ( kind = dp ) :: mySigma_i,mySigma_j
156 <    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
155 >    real ( kind = dp ) :: Sigma_i, Sigma_j
156 >    real ( kind = dp ) :: Epsilon_i, Epsilon_j
157      real ( kind = dp ) :: rcut6
158 <    logical :: I_isLJ, J_isLJ
158 >
159      status = 0
160      
161 <    nAtypes = getSize(atypes)
161 >    nAtypes = size(ParameterMap)
162 >    
163      if (nAtypes == 0) then
164         status = -1
165         return
166      end if
133        
134    if (.not. associated(ljMixed)) then
135       allocate(ljMixed(nAtypes, nAtypes))
136    endif
167  
168 <    rcut6 = LJ_rcut**6
169 <
170 < ! This loops through all atypes, even those that don't support LJ forces.
168 >    if (.not.have_rcut) then
169 >       status = -1
170 >       return
171 >    endif
172 >    
173 >    if (.not. allocated(MixingMap)) then
174 >       allocate(MixingMap(nAtypes, nAtypes))
175 >    endif
176 >    
177 >    rcut6 = LJ_rcut**6
178 >    
179 >    ! This loops through all atypes, even those that don't support LJ forces.
180      do i = 1, nAtypes
181 <
182 <       call getElementProperty(atypes, i, "is_LJ", I_isLJ)
183 <
184 <       if (I_isLJ) then
185 <
186 <          call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
187 <          call getElementProperty(atypes, i, "lj_sigma",   mySigma_i)
188 <          ! do self mixing rule
189 <          ljMixed(i,i)%sigma   = mySigma_i
181 >      
182 >       Epsilon_i = ParameterMap(i)%epsilon
183 >       Sigma_i = ParameterMap(i)%sigma
184 >      
185 >       ! do self mixing rule
186 >       MixingMap(i,i)%sigma   = Sigma_i          
187 >       MixingMap(i,i)%sigma6  = Sigma_i ** 6          
188 >       MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
189 >       MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
190 >       MixingMap(i,i)%epsilon = Epsilon_i          
191 >       MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
192 >            (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
193 >      
194 >       do j = i + 1, nAtypes
195            
196 <          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
196 >          Epsilon_j = ParameterMap(j)%epsilon
197 >          Sigma_j = ParameterMap(j)%sigma
198            
199 <          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
199 >          ! only the distance parameter uses different mixing policies
200 >          if (useGeometricDistanceMixing) then
201 >             ! only for OPLS as far as we can tell
202 >             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
203 >          else
204 >             ! everyone else
205 >             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
206 >          endif
207            
208 <          ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
208 >          ! energy parameter is always geometric mean:
209 >          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
210 >                    
211 >          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
212 >          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
213 >          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
214            
215 +          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
216 +               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
217            
218 <          ljMixed(i,i)%epsilon = myEpsilon_i
218 >          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
219 >          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
220 >          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
221 >          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
222 >          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
223 >          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
224            
225 <          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
225 >       end do
226      end do
227      
228 <  end subroutine createMixingList
229 <  
228 >  end subroutine createMixingMap
229 >        
230    subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
231         pot, f, do_pot)
232  
# Line 225 | Line 248 | contains
248      real( kind = dp ) :: t6
249      real( kind = dp ) :: t12
250      real( kind = dp ) :: delta
251 <    integer :: id1, id2
251 >    integer :: id1, id2, localError
252  
253 +    if (.not.haveMixingMap) then
254 +       localError = 0
255 +       call createMixingMap(localError)
256 +       if ( localError .ne. 0 ) then
257 +          call handleError("LJ", "MixingMap creation failed!")
258 +          return
259 +       end if
260 +    endif
261 +
262      ! Look up the correct parameters in the mixing matrix
263   #ifdef IS_MPI
264 <    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
265 <    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
266 <    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
264 >    sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
265 >    epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
266 >    delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
267   #else
268 <    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
269 <    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
270 <    delta    = ljMixed(atid(atom1),atid(atom2))%delta
268 >    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
269 >    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
270 >    delta    = MixingMap(atid(atom1),atid(atom2))%delta
271   #endif
272  
273      r6 = r2 * r2 * r2
# Line 309 | Line 341 | contains
341    
342    
343    !! Calculates the mixing for sigma or epslon
344 +    
345 + end module lj
346 +
347 + subroutine newLJtype(ident, sigma, epsilon, status)
348 +  use lj, ONLY : module_newLJtype => newLJtype
349 +  integer, parameter :: DP = selected_real_kind(15)
350 +  integer,intent(inout) :: ident
351 +  real(kind=dp),intent(inout) :: sigma
352 +  real(kind=dp),intent(inout) :: epsilon
353 +  integer,intent(inout) :: status
354    
355 <  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
356 <    character(len=*) :: thisParam
357 <    real(kind = dp)  :: param1
316 <    real(kind = dp)  :: param2
317 <    real(kind = dp ) :: myMixParam
355 >  call module_newLJtype(ident, sigma, epsilon, status)
356 >  
357 > end subroutine newLJtype
358  
359 <    integer, optional :: status  
360 <
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
359 > subroutine useGeometricMixing()
360 >  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
361    
362 < end module lj
362 >  call module_useGeometricMixing()
363 >  return
364 > end subroutine useGeometricMixing

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines