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 2717 by gezelter, Mon Apr 17 21:49:12 2006 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.21 2006-04-17 21:49:12 gezelter Exp $, $Date: 2006-04-17 21:49:12 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
47  
48 +
49   module lj
8 –  use definitions
50    use atype_module
10 –  use switcheroo
51    use vector_class
52    use simulation
53 +  use status
54 +  use fForceOptions
55 +  use interpolation
56   #ifdef IS_MPI
57    use mpiSimulation
58   #endif
# Line 17 | Line 60 | module lj
60  
61    implicit none
62    PRIVATE
20 –
63   #define __FORTRAN90
64 < #include "UseTheForce/fForceField.h"
64 > #include "UseTheForce/DarkSide/fInteractionMap.h"
65  
66 <  integer, save :: LJ_Mixing_Policy
67 <  real(kind=DP), save :: LJ_rcut
68 <  logical, save :: havePolicy = .false.
69 <  logical, save :: haveCut = .false.
70 <  logical, save :: LJ_do_shift = .false.
71 <  
72 <  !! Logical has lj force field module been initialized?
73 <  
74 <  logical, save :: LJ_FF_initialized = .false.
75 <  
76 <  !! Public methods and data
77 <  public :: init_LJ_FF
78 <  public :: setCutoffLJ
66 >  integer, parameter :: DP = selected_real_kind(15)
67 >
68 >  logical, save :: useGeometricDistanceMixing = .false.
69 >  logical, save :: haveMixingMap = .false.
70 >  logical, save :: useSplines = .false.
71 >
72 >  real(kind=DP), save :: defaultCutoff = 0.0_DP
73 >  logical, save :: defaultShift = .false.
74 >  logical, save :: haveDefaultCutoff = .false.
75 >
76 >  type, private :: LJtype
77 >     integer       :: atid
78 >     real(kind=dp) :: sigma
79 >     real(kind=dp) :: epsilon
80 >     logical       :: isSoftCore = .false.
81 >  end type LJtype
82 >
83 >  type, private :: LJList
84 >     integer               :: Nljtypes = 0
85 >     integer               :: currentLJtype = 0
86 >     type(LJtype), pointer :: LJtypes(:)      => null()
87 >     integer, pointer      :: atidToLJtype(:) => null()
88 >  end type LJList
89 >
90 >  type(LJList), save :: LJMap
91 >
92 >  type :: MixParameters
93 >     real(kind=DP) :: sigma
94 >     real(kind=DP) :: epsilon
95 >     real(kind=dp) :: sigmai
96 >     real(kind=dp) :: rCut
97 >     logical       :: rCutWasSet = .false.
98 >     logical       :: shiftedPot
99 >     logical       :: isSoftCore = .false.
100 >  end type MixParameters
101 >
102 >  type(MixParameters), dimension(:,:), allocatable :: MixingMap
103 >
104 >  type(cubicSpline), save :: vLJspline
105 >  type(cubicSpline), save :: vLJpspline
106 >  type(cubicSpline), save :: vSoftSpline
107 >  type(cubicSpline), save :: vSoftpSpline
108 >
109 >  public :: newLJtype
110 >  public :: setLJDefaultCutoff
111 >  public :: getSigma
112 >  public :: getEpsilon
113    public :: do_lj_pair
114 <  
115 <  type :: lj_mixed_params
116 <     !! 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 <  
114 >  public :: destroyLJtypes
115 >  public :: setLJsplineRmax
116 >
117   contains
118  
119 <  subroutine init_LJ_FF(mix_Policy, status)
120 <    integer, intent(in) :: mix_Policy
121 <    integer, intent(out) :: status
122 <    integer :: myStatus
123 <    
124 <    if (mix_Policy == LB_MIXING_RULE) then
125 <       LJ_Mixing_Policy = LB_MIXING_RULE
119 >  subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
120 >    integer,intent(in) :: c_ident
121 >    real(kind=dp),intent(in) :: sigma
122 >    real(kind=dp),intent(in) :: epsilon
123 >    integer, intent(in) :: isSoftCore
124 >    integer,intent(out) :: status
125 >    integer :: nLJTypes, ntypes, myATID
126 >    integer, pointer :: MatchList(:) => null()
127 >    integer :: current
128 >
129 >    status = 0
130 >    ! check to see if this is the first time into this routine...
131 >    if (.not.associated(LJMap%LJtypes)) then
132 >
133 >       call getMatchingElementList(atypes, "is_LennardJones", .true., &
134 >            nLJTypes, MatchList)
135 >      
136 >       LJMap%nLJtypes =  nLJTypes
137 >
138 >       allocate(LJMap%LJtypes(nLJTypes))
139 >
140 >       ntypes = getSize(atypes)
141 >
142 >       allocate(LJMap%atidToLJtype(ntypes))
143 >    end if
144 >
145 >    LJMap%currentLJtype = LJMap%currentLJtype + 1
146 >    current = LJMap%currentLJtype
147 >
148 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
149 >    LJMap%atidToLJtype(myATID)        = current
150 >    LJMap%LJtypes(current)%atid       = myATID
151 >    LJMap%LJtypes(current)%sigma      = sigma
152 >    LJMap%LJtypes(current)%epsilon    = epsilon
153 >    if (isSoftCore .eq. 1) then
154 >       LJMap%LJtypes(current)%isSoftCore = .true.
155      else
156 <       if (mix_Policy == EXPLICIT_MIXING_RULE) then
65 <          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
66 <       else
67 <          write(*,*) 'Unknown Mixing Policy!'
68 <          status = -1
69 <          return
70 <       endif
156 >       LJMap%LJtypes(current)%isSoftCore = .false.
157      endif
158 +  end subroutine newLJtype
159  
160 <    havePolicy = .true.
160 >  subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
161 >    real(kind=dp), intent(in) :: thisRcut
162 >    logical, intent(in) :: shiftedPot
163 >    defaultCutoff = thisRcut
164 >    defaultShift = shiftedPot
165 >    haveDefaultCutoff = .true.
166 >    !we only want to build LJ Mixing map and spline if LJ is being used.
167 >    if(LJMap%nLJTypes /= 0) then
168 >       call createMixingMap()
169 >       call setLJsplineRmax(defaultCutoff)
170 >    end if
171  
172 <    if (haveCut) then
173 <       status = 0
174 <       call createMixingList(myStatus)
175 <       if (myStatus /= 0) then
176 <          status = -1
177 <          return
178 <       end if
179 <      
180 <       LJ_FF_initialized = .true.
172 >  end subroutine setLJDefaultCutoff
173 >
174 >  function getSigma(atid) result (s)
175 >    integer, intent(in) :: atid
176 >    integer :: ljt1
177 >    real(kind=dp) :: s
178 >
179 >    if (LJMap%currentLJtype == 0) then
180 >       call handleError("LJ", "No members in LJMap")
181 >       return
182      end if
85 –  
86 –  end subroutine init_LJ_FF
87 –  
88 –  subroutine setCutoffLJ(rcut, do_shift, status)
89 –    logical, intent(in):: do_shift
90 –    integer :: status, myStatus
91 –    real(kind=dp) :: rcut
183  
184 < #define __FORTRAN90
185 < #include "UseTheForce/fSwitchingFunction.h"
184 >    ljt1 = LJMap%atidToLJtype(atid)
185 >    s = LJMap%LJtypes(ljt1)%sigma
186  
187 <    status = 0
187 >  end function getSigma
188  
189 <    LJ_rcut = rcut
190 <    LJ_do_shift = do_shift
191 <    call set_switch(LJ_SWITCH, rcut, rcut)
192 <    haveCut = .true.
189 >  function getEpsilon(atid) result (e)
190 >    integer, intent(in) :: atid
191 >    integer :: ljt1
192 >    real(kind=dp) :: e
193  
194 <    if (havePolicy) then
195 <       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    
113 <    
114 <    return
115 <  end subroutine setCutoffLJ
116 <  
117 <  subroutine createMixingList(status)
118 <    integer :: nAtypes
119 <    integer :: status
120 <    integer :: i
121 <    integer :: j
122 <    real ( kind = dp ) :: mySigma_i,mySigma_j
123 <    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
124 <    real ( kind = dp ) :: rcut6
125 <    logical :: I_isLJ, J_isLJ
126 <    status = 0
127 <    
128 <    nAtypes = getSize(atypes)
129 <    if (nAtypes == 0) then
130 <       status = -1
194 >    if (LJMap%currentLJtype == 0) then
195 >       call handleError("LJ", "No members in LJMap")
196         return
197      end if
133 –        
134 –    if (.not. associated(ljMixed)) then
135 –       allocate(ljMixed(nAtypes, nAtypes))
136 –    endif
198  
199 <    rcut6 = LJ_rcut**6
199 >    ljt1 = LJMap%atidToLJtype(atid)
200 >    e = LJMap%LJtypes(ljt1)%epsilon
201  
202 < ! This loops through all atypes, even those that don't support LJ forces.
141 <    do i = 1, nAtypes
202 >  end function getEpsilon
203  
204 <       call getElementProperty(atypes, i, "is_LJ", I_isLJ)
204 >  subroutine createMixingMap()
205 >    integer :: nLJtypes, i, j
206 >    real ( kind = dp ) :: s1, s2, e1, e2
207 >    real ( kind = dp ) :: rcut6, tp6, tp12
208 >    logical :: isSoftCore1, isSoftCore2, doShift
209  
210 <       if (I_isLJ) then
210 >    if (LJMap%currentLJtype == 0) then
211 >       call handleError("LJ", "No members in LJMap")
212 >       return
213 >    end if
214  
215 <          call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
216 <          call getElementProperty(atypes, i, "lj_sigma",   mySigma_i)
217 <          ! do self mixing rule
218 <          ljMixed(i,i)%sigma   = mySigma_i
215 >    nLJtypes = LJMap%nLJtypes
216 >
217 >    if (.not. allocated(MixingMap)) then
218 >       allocate(MixingMap(nLJtypes, nLJtypes))
219 >    endif
220 >
221 >    useGeometricDistanceMixing = usesGeometricDistanceMixing()
222 >    do i = 1, nLJtypes
223 >
224 >       s1 = LJMap%LJtypes(i)%sigma
225 >       e1 = LJMap%LJtypes(i)%epsilon
226 >       isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
227 >
228 >       do j = i, nLJtypes
229            
230 <          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
230 >          s2 = LJMap%LJtypes(j)%sigma
231 >          e2 = LJMap%LJtypes(j)%epsilon
232 >          isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
233            
234 <          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
155 <          
156 <          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)
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
203 <    end do
204 <    
205 <  end subroutine createMixingList
206 <  
207 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
208 <       pot, f, do_pot)
234 >          MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
235  
236 +          ! only the distance parameter uses different mixing policies
237 +          if (useGeometricDistanceMixing) then
238 +             MixingMap(i,j)%sigma = dsqrt(s1 * s2)
239 +          else
240 +             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
241 +          endif
242 +          
243 +          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
244 +
245 +          MixingMap(i,j)%sigmai = 1.0_DP  / (MixingMap(i,j)%sigma)
246 +
247 +          if (haveDefaultCutoff) then
248 +             MixingMap(i,j)%shiftedPot = defaultShift
249 +          else
250 +             MixingMap(i,j)%shiftedPot = defaultShift
251 +          endif          
252 +
253 +          if (i.ne.j) then
254 +             MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
255 +             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
256 +             MixingMap(j,i)%sigmai     = MixingMap(i,j)%sigmai
257 +             MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
258 +             MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
259 +             MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
260 +             MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
261 +          endif
262 +
263 +       enddo
264 +    enddo
265 +    
266 +    haveMixingMap = .true.
267 +    
268 +  end subroutine createMixingMap
269 +  
270 +  subroutine setLJsplineRmax(largestRcut)
271 +    real( kind = dp ), intent(in) :: largestRcut
272 +    real( kind = dp ) :: s, bigS, smallS, rmax, rmin
273 +    integer :: np, i
274 +
275 +    if (LJMap%nLJtypes .ne. 0) then
276 +
277 +       !
278 +       ! find the largest and smallest values of sigma that we'll need
279 +       !
280 +       bigS = 0.0_DP
281 +       smallS = 1.0e9
282 +       do i = 1, LJMap%nLJtypes
283 +          s = LJMap%LJtypes(i)%sigma
284 +          if (s .gt. bigS) bigS = s
285 +          if (s .lt. smallS) smallS = s
286 +       enddo
287 +      
288 +       !
289 +       ! give ourselves a 20% margin just in case
290 +       !
291 +       rmax = 1.2 * largestRcut / smallS    
292 +       !
293 +       ! assume atoms will never get closer than 1 angstrom
294 +       !
295 +       rmin = 1 / bigS
296 +       !
297 +       ! assume 500 points is enough
298 +       !
299 +       np = 500
300 +      
301 +       write(*,*) 'calling setupSplines with rmin = ', rmin, ' rmax = ', rmax, &
302 +            ' np = ', np
303 +      
304 +       call setupSplines(rmin, rmax, np)
305 +
306 +    endif
307 +    return
308 +  end subroutine setLJsplineRmax
309 +        
310 +  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
311 +       pot, f, do_pot)
312 +    
313      integer, intent(in) ::  atom1, atom2
314 <    real( kind = dp ), intent(in) :: rij, r2
314 >    integer :: atid1, atid2, ljt1, ljt2
315 >    real( kind = dp ), intent(in) :: rij, r2, rcut
316      real( kind = dp ) :: pot, sw, vpair
317      real( kind = dp ), dimension(3,nLocal) :: f    
318      real( kind = dp ), intent(in), dimension(3) :: d
# Line 218 | Line 322 | contains
322      ! local Variables
323      real( kind = dp ) :: drdx, drdy, drdz
324      real( kind = dp ) :: fx, fy, fz
325 +    real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
326      real( kind = dp ) :: pot_temp, dudr
327 <    real( kind = dp ) :: sigma6
327 >    real( kind = dp ) :: sigmai
328      real( kind = dp ) :: epsilon
329 <    real( kind = dp ) :: r6
330 <    real( kind = dp ) :: t6
226 <    real( kind = dp ) :: t12
227 <    real( kind = dp ) :: delta
228 <    integer :: id1, id2
329 >    logical :: isSoftCore, shiftedPot
330 >    integer :: id1, id2, localError
331  
332 +    if (.not.haveMixingMap) then
333 +       call createMixingMap()
334 +    endif
335 +
336      ! Look up the correct parameters in the mixing matrix
337   #ifdef IS_MPI
338 <    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
339 <    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
234 <    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
338 >    atid1 = atid_Row(atom1)
339 >    atid2 = atid_Col(atom2)
340   #else
341 <    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
342 <    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
238 <    delta    = ljMixed(atid(atom1),atid(atom2))%delta
341 >    atid1 = atid(atom1)
342 >    atid2 = atid(atom2)
343   #endif
344  
345 <    r6 = r2 * r2 * r2
346 <    
347 <    t6  = sigma6/ r6
348 <    t12 = t6 * t6    
349 <  
350 <    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
351 <    if (LJ_do_shift) then
352 <       pot_temp = pot_temp + delta
345 >    ljt1 = LJMap%atidToLJtype(atid1)
346 >    ljt2 = LJMap%atidToLJtype(atid2)
347 >
348 >    sigmai     = MixingMap(ljt1,ljt2)%sigmai
349 >    epsilon    = MixingMap(ljt1,ljt2)%epsilon
350 >    isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
351 >    shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
352 >
353 >    ros = rij * sigmai
354 >    myPotC = 0.0_DP
355 >
356 >    if (isSoftCore) then
357 >
358 >       call getSoftFunc(ros, myPot, myDeriv)
359 >
360 >       if (shiftedPot) then
361 >          rcos = rcut * sigmai
362 >          call getSoftFunc(rcos, myPotC, myDerivC)
363 >       endif
364 >              
365 >    else
366 >
367 >       call getLJfunc(ros, myPot, myDeriv)
368 >
369 >       if (shiftedPot) then
370 >          rcos = rcut * sigmai
371 >          call getLJfunc(rcos, myPotC, myDerivC)
372 >       endif
373 >      
374      endif
375  
376 +    !write(*,*) rij, ros, rcos, myPot, myDeriv, myPotC
377 +
378 +    pot_temp = epsilon * (myPot - myPotC)
379      vpair = vpair + pot_temp
380 <      
381 <    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
254 <      
380 >    dudr = sw * epsilon * myDeriv * sigmai
381 >
382      drdx = d(1) / rij
383      drdy = d(2) / rij
384      drdz = d(3) / rij
385 <      
385 >
386      fx = dudr * drdx
387      fy = dudr * drdy
388      fz = dudr * drdz
389 <    
263 <      
389 >
390   #ifdef IS_MPI
391      if (do_pot) then
392 <       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
393 <       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
392 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
393 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
394      endif
395 <    
395 >
396      f_Row(1,atom1) = f_Row(1,atom1) + fx
397      f_Row(2,atom1) = f_Row(2,atom1) + fy
398      f_Row(3,atom1) = f_Row(3,atom1) + fz
399 <    
399 >
400      f_Col(1,atom2) = f_Col(1,atom2) - fx
401      f_Col(2,atom2) = f_Col(2,atom2) - fy
402      f_Col(3,atom2) = f_Col(3,atom2) - fz      
403 <    
403 >
404   #else
405      if (do_pot) pot = pot + sw*pot_temp
406  
407      f(1,atom1) = f(1,atom1) + fx
408      f(2,atom1) = f(2,atom1) + fy
409      f(3,atom1) = f(3,atom1) + fz
410 <    
410 >
411      f(1,atom2) = f(1,atom2) - fx
412      f(2,atom2) = f(2,atom2) - fy
413      f(3,atom2) = f(3,atom2) - fz
414   #endif
415 <        
415 >
416   #ifdef IS_MPI
417      id1 = AtomRowToGlobal(atom1)
418      id2 = AtomColToGlobal(atom2)
# Line 296 | Line 422 | contains
422   #endif
423  
424      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
425 <      
425 >
426         fpair(1) = fpair(1) + fx
427         fpair(2) = fpair(2) + fy
428         fpair(3) = fpair(3) + fz
# Line 304 | Line 430 | contains
430      endif
431  
432      return    
433 <    
433 >
434    end subroutine do_lj_pair
309 –  
310 –  
311 –  !! 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
435  
436 <    integer, optional :: status  
436 >  subroutine destroyLJTypes()
437  
438 <    myMixParam = 0.0_dp
438 >    LJMap%nLJtypes = 0
439 >    LJMap%currentLJtype = 0
440      
441 <    if (present(status)) status = 0
442 <    select case (LJ_Mixing_Policy)
443 <    case (1)
444 <       select case (thisParam)
445 <       case ("sigma")
446 <          myMixParam = 0.5_dp * (param1 + param2)
447 <       case ("epsilon")
448 <          myMixParam = sqrt(param1 * param2)
449 <       case default
450 <          status = -1
451 <       end select
452 <    case default
453 <       status = -1
454 <    end select
455 <  end function calcLJMix
456 <  
441 >    if (associated(LJMap%LJtypes)) then
442 >       deallocate(LJMap%LJtypes)
443 >       LJMap%LJtypes => null()
444 >    end if
445 >    
446 >    if (associated(LJMap%atidToLJtype)) then
447 >       deallocate(LJMap%atidToLJtype)
448 >       LJMap%atidToLJtype => null()
449 >    end if
450 >    
451 >    haveMixingMap = .false.
452 >
453 >    call deleteSpline(vLJspline)
454 >    call deleteSpline(vLJpspline)
455 >    call deleteSpline(vSoftSpline)
456 >    call deleteSpline(vSoftpSpline)
457 >
458 >  end subroutine destroyLJTypes
459 >
460 >  subroutine getLJfunc(r, myPot, myDeriv)
461 >
462 >    real(kind=dp), intent(in) :: r
463 >    real(kind=dp), intent(inout) :: myPot, myDeriv
464 >    real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
465 >    real(kind=dp) :: a, b, c, d, dx
466 >    integer :: j
467 >
468 >    if (useSplines) then
469 >       j = MAX(1, MIN(vLJspline%np, idint((r-vLJspline%x(1)) * vLJspline%dx_i) + 1))
470 >      
471 >       dx = r - vLJspline%x(j)
472 >      
473 >       a = vLJspline%c(1,j)
474 >       b = vLJspline%c(2,j)
475 >       c = vLJspline%c(3,j)
476 >       d = vLJspline%c(4,j)
477 >      
478 >       myPot = c + dx * d
479 >       myPot = b + dx * myPot
480 >       myPot = a + dx * myPot
481 >
482 >       a = vLJpspline%c(1,j)
483 >       b = vLJpspline%c(2,j)
484 >       c = vLJpspline%c(3,j)
485 >       d = vLJpspline%c(4,j)
486 >      
487 >       myDeriv = c + dx * d
488 >       myDeriv = b + dx * myDeriv  
489 >       myDeriv = a + dx * myDeriv
490 >      
491 >    else
492 >       ri = 1.0_DP / r
493 >       ri2 = ri*ri
494 >       ri6 = ri2*ri2*ri2
495 >       ri7 = ri6*ri
496 >       ri12 = ri6*ri6
497 >       ri13 = ri12*ri
498 >      
499 >       myPot = 4.0_DP * (ri12 - ri6)
500 >       myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
501 >    endif
502 >
503 >    return
504 >  end subroutine getLJfunc
505 >
506 >  subroutine getSoftFunc(r, myPot, myDeriv)
507 >    
508 >    real(kind=dp), intent(in) :: r
509 >    real(kind=dp), intent(inout) :: myPot, myDeriv
510 >    real(kind=dp) :: ri, ri2, ri6, ri7
511 >    real(kind=dp) :: a, b, c, d, dx
512 >    integer :: j
513 >    
514 >    if (useSplines) then
515 >       j = MAX(1, MIN(vSoftSpline%np, idint((r-vSoftSpline%x(1)) * vSoftSpline%dx_i) + 1))
516 >      
517 >       dx = r - vSoftSpline%x(j)
518 >      
519 >       a = vSoftSpline%c(1,j)
520 >       b = vSoftSpline%c(2,j)
521 >       c = vSoftSpline%c(3,j)
522 >       d = vSoftSpline%c(4,j)
523 >      
524 >       myPot = c + dx * d
525 >       myPot = b + dx * myPot
526 >       myPot = a + dx * myPot
527 >
528 >       a = vSoftPspline%c(1,j)
529 >       b = vSoftPspline%c(2,j)
530 >       c = vSoftPspline%c(3,j)
531 >       d = vSoftPspline%c(4,j)
532 >      
533 >       myDeriv = c + dx * d
534 >       myDeriv = b + dx * myDeriv  
535 >       myDeriv = a + dx * myDeriv
536 >      
537 >    else
538 >       ri = 1.0_DP / r    
539 >       ri2 = ri*ri
540 >       ri6 = ri2*ri2*ri2
541 >       ri7 = ri6*ri
542 >       myPot = 4.0_DP * (ri6)
543 >       myDeriv = - 24.0_DP * ri7
544 >    endif
545 >
546 >    return
547 >  end subroutine getSoftFunc
548 >
549 >  subroutine setupSplines(rmin, rmax, np)
550 >    real( kind = dp ), intent(in) :: rmin, rmax
551 >    integer, intent(in) :: np
552 >    real( kind = dp ) :: rvals(np), vLJ(np), vLJp(np), vSoft(np), vSoftp(np)
553 >    real( kind = dp ) :: dr, r, ri, ri2, ri6, ri7, ri12, ri13
554 >    real( kind = dp ) :: vljpp1, vljppn, vsoftpp1, vsoftppn
555 >    integer :: i
556 >
557 >    dr = (rmax-rmin) / float(np-1)
558 >    
559 >    do i = 1, np
560 >       r = rmin + dble(i-1)*dr
561 >       ri = 1.0_DP / r
562 >       ri2 = ri*ri
563 >       ri6 = ri2*ri2*ri2
564 >       ri7 = ri6*ri
565 >       ri12 = ri6*ri6
566 >       ri13 = ri12*ri
567 >
568 >       rvals(i) = r
569 >       vLJ(i) = 4.0_DP * (ri12 - ri6)
570 >       vLJp(i) = 24.0_DP * (ri7 - 2.0_DP * ri13)
571 >
572 >       vSoft(i) = 4.0_DP * (ri6)
573 >       vSoftp(i) = - 24.0_DP * ri7      
574 >    enddo
575 >
576 >    vljpp1 = 624.0_DP / (rmin)**(14)  - 168.0_DP / (rmin)**(8)
577 >    vljppn = 624.0_DP / (rmax)**(14)  - 168.0_DP / (rmax)**(8)
578 >
579 >    vsoftpp1 = 168.0_DP / (rmin)**(8)
580 >    vsoftppn = 168.0_DP / (rmax)**(8)
581 >
582 >    call newSpline(vLJspline, rvals, vLJ, vLJp(1), vLJp(np), .true.)
583 >    call newSpline(vLJpspline, rvals, vLJp, vljpp1, vljppn, .true.)
584 >    call newSpline(vSoftSpline, rvals, vSoft, vSoftp(1), vSoftp(np), .true.)
585 >    call newSpline(vSoftpSpline, rvals, vSoftp, vsoftpp1, vsoftppn, .true.)
586 >
587 >    return
588 >  end subroutine setupSplines
589   end module lj

Diff Legend

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