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 1633 by gezelter, Fri Oct 22 20:22:48 2004 UTC vs.
Revision 2722 by gezelter, Thu Apr 20 18:24:24 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.
44   !! @author Charles F. Vardeman II
45   !! @author Matthew Meineke
46 < !! @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 $
46 > !! @version $Id: LJ.F90,v 1.22 2006-04-20 18:24:24 gezelter Exp $, $Date: 2006-04-20 18:24:24 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
47  
48 +
49   module lj
50    use atype_module
8 –  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 16 | Line 60 | module lj
60  
61    implicit none
62    PRIVATE
63 <  
63 > #define __FORTRAN90
64 > #include "UseTheForce/DarkSide/fInteractionMap.h"
65 >
66    integer, parameter :: DP = selected_real_kind(15)
67 <  
68 <  type, private :: LjType
69 <     integer :: ident
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 <  end type LjType
81 <  
82 <  type(LjType), dimension(:), allocatable :: ParameterMap
83 <  
84 <  logical, save :: haveMixingMap = .false.
85 <  
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)  :: sigma6
96 <     real(kind=dp)  :: tp6
97 <     real(kind=dp)  :: tp12
98 <     real(kind=dp)  :: delta
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 <  
101 >
102    type(MixParameters), dimension(:,:), allocatable :: MixingMap
103 <  
104 <  real(kind=DP), save :: LJ_rcut
105 <  logical, save :: have_rcut = .false.
106 <  logical, save :: LJ_do_shift = .false.
107 <  logical, save :: useGeometricDistanceMixing = .false.
108 <  
109 <  !! Public methods and data
110 <  
50 <  public :: setCutoffLJ
51 <  public :: useGeometricMixing
52 <  public :: do_lj_pair
53 <  public :: newLJtype  
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 <  
113 >  public :: do_lj_pair
114 >  public :: destroyLJtypes
115 >  public :: setLJsplineRmax
116 >
117   contains
118  
119 <  subroutine newLJtype(ident, sigma, epsilon, status)
120 <    integer,intent(in) :: ident
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 :: nAtypes
125 >    integer :: nLJTypes, ntypes, myATID
126 >    integer, pointer :: MatchList(:) => null()
127 >    integer :: current
128  
129      status = 0
130 <    
131 <    !! Be simple-minded and assume that we need a ParameterMap that
69 <    !! is the same size as the total number of atom types
130 >    ! check to see if this is the first time into this routine...
131 >    if (.not.associated(LJMap%LJtypes)) then
132  
133 <    if (.not.allocated(ParameterMap)) then
133 >       call getMatchingElementList(atypes, "is_LennardJones", .true., &
134 >            nLJTypes, MatchList)
135        
136 <       nAtypes = getSize(atypes)
137 <    
138 <       if (nAtypes == 0) then
139 <          status = -1
140 <          return
141 <       end if
142 <      
80 <       if (.not. allocated(ParameterMap)) then
81 <          allocate(ParameterMap(nAtypes))
82 <       endif
83 <      
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 <    if (ident .gt. size(ParameterMap)) then
146 <       status = -1
88 <       return
89 <    endif
90 <    
91 <    ! set the values for ParameterMap for this atom type:
145 >    LJMap%currentLJtype = LJMap%currentLJtype + 1
146 >    current = LJMap%currentLJtype
147  
148 <    ParameterMap(ident)%ident = ident
149 <    ParameterMap(ident)%epsilon = epsilon
150 <    ParameterMap(ident)%sigma = sigma
151 <    
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 >       LJMap%LJtypes(current)%isSoftCore = .false.
157 >    endif
158    end subroutine newLJtype
159  
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 +  end subroutine setLJDefaultCutoff
173 +
174    function getSigma(atid) result (s)
175      integer, intent(in) :: atid
176 <    integer :: localError
176 >    integer :: ljt1
177      real(kind=dp) :: s
178 <    
179 <    if (.not.allocated(ParameterMap)) then
180 <       call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
178 >
179 >    if (LJMap%currentLJtype == 0) then
180 >       call handleError("LJ", "No members in LJMap")
181         return
182      end if
183 <    
184 <    s = ParameterMap(atid)%sigma
183 >
184 >    ljt1 = LJMap%atidToLJtype(atid)
185 >    s = LJMap%LJtypes(ljt1)%sigma
186 >
187    end function getSigma
188  
189    function getEpsilon(atid) result (e)
190      integer, intent(in) :: atid
191 <    integer :: localError
191 >    integer :: ljt1
192      real(kind=dp) :: e
193 <    
194 <    if (.not.allocated(ParameterMap)) then
195 <       call handleError("dipole-dipole", "no ParameterMap was present before first call of getEpsilon!")
193 >
194 >    if (LJMap%currentLJtype == 0) then
195 >       call handleError("LJ", "No members in LJMap")
196         return
197      end if
121 –    
122 –    e = ParameterMap(atid)%epsilon
123 –  end function getEpsilon
198  
199 +    ljt1 = LJMap%atidToLJtype(atid)
200 +    e = LJMap%LJtypes(ljt1)%epsilon
201  
202 <  subroutine setCutoffLJ(rcut, do_shift, status)
127 <    logical, intent(in):: do_shift
128 <    integer :: status, myStatus
129 <    real(kind=dp) :: rcut
202 >  end function getEpsilon
203  
204 < #define __FORTRAN90
205 < #include "UseTheForce/fSwitchingFunction.h"
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 <    status = 0
211 <
136 <    LJ_rcut = rcut
137 <    LJ_do_shift = do_shift
138 <    call set_switch(LJ_SWITCH, rcut, rcut)
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 createMixingMap(status)
151 <    integer :: nAtypes
152 <    integer :: status
153 <    integer :: i
154 <    integer :: j
155 <    real ( kind = dp ) :: Sigma_i, Sigma_j
156 <    real ( kind = dp ) :: Epsilon_i, Epsilon_j
157 <    real ( kind = dp ) :: rcut6
158 <
159 <    status = 0
160 <    
161 <    nAtypes = size(ParameterMap)
162 <    
163 <    if (nAtypes == 0) then
164 <       status = -1
210 >    if (LJMap%currentLJtype == 0) then
211 >       call handleError("LJ", "No members in LJMap")
212         return
213      end if
214  
215 <    if (.not.have_rcut) then
216 <       status = -1
170 <       return
171 <    endif
172 <    
215 >    nLJtypes = LJMap%nLJtypes
216 >
217      if (.not. allocated(MixingMap)) then
218 <       allocate(MixingMap(nAtypes, nAtypes))
218 >       allocate(MixingMap(nLJtypes, nLJtypes))
219      endif
220 <    
221 <    rcut6 = LJ_rcut**6
222 <    
223 <    ! This loops through all atypes, even those that don't support LJ forces.
224 <    do i = 1, nAtypes
225 <      
226 <       Epsilon_i = ParameterMap(i)%epsilon
227 <       Sigma_i = ParameterMap(i)%sigma
228 <      
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
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 <          Epsilon_j = ParameterMap(j)%epsilon
231 <          Sigma_j = ParameterMap(j)%sigma
230 >          s2 = LJMap%LJtypes(j)%sigma
231 >          e2 = LJMap%LJtypes(j)%epsilon
232 >          isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
233            
234 +          MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
235 +
236            ! only the distance parameter uses different mixing policies
237            if (useGeometricDistanceMixing) then
238 <             ! only for OPLS as far as we can tell
202 <             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
238 >             MixingMap(i,j)%sigma = dsqrt(s1 * s2)
239            else
240 <             ! everyone else
205 <             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
240 >             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
241            endif
242            
243 <          ! energy parameter is always geometric mean:
244 <          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
245 <                    
246 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
247 <          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
248 <          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
249 <          
250 <          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
251 <               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
252 <          
253 <          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
254 <          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
255 <          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
256 <          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
257 <          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
258 <          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
259 <          
260 <       end do
261 <    end do
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, sw, vpair, fpair, &
310 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
311         pot, f, do_pot)
312 <
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 241 | 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
248 <    real( kind = dp ) :: t6
249 <    real( kind = dp ) :: t12
250 <    real( kind = dp ) :: delta
329 >    logical :: isSoftCore, shiftedPot
330      integer :: id1, id2, localError
331  
332      if (.not.haveMixingMap) then
333 <       localError = 0
255 <       call createMixingMap(localError)
256 <       if ( localError .ne. 0 ) then
257 <          call handleError("LJ", "MixingMap creation failed!")
258 <          return
259 <       end if
333 >       call createMixingMap()
334      endif
335  
336      ! Look up the correct parameters in the mixing matrix
337   #ifdef IS_MPI
338 <    sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
339 <    epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
266 <    delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
338 >    atid1 = atid_Row(atom1)
339 >    atid2 = atid_Col(atom2)
340   #else
341 <    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
342 <    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
270 <    delta    = MixingMap(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
286 <      
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 <    
295 <      
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 328 | 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 336 | Line 430 | contains
430      endif
431  
432      return    
433 <    
433 >
434    end subroutine do_lj_pair
435 <  
436 <  
437 <  !! Calculates the mixing for sigma or epslon
435 >
436 >  subroutine destroyLJTypes()
437 >
438 >    LJMap%nLJtypes = 0
439 >    LJMap%currentLJtype = 0
440      
441 < end module lj
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 < subroutine newLJtype(ident, sigma, epsilon, status)
454 <  use lj, ONLY : module_newLJtype => newLJtype
455 <  integer, parameter :: DP = selected_real_kind(15)
456 <  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 <  call module_newLJtype(ident, sigma, epsilon, status)
356 <  
357 < end subroutine newLJtype
453 >    call deleteSpline(vLJspline)
454 >    call deleteSpline(vLJpspline)
455 >    call deleteSpline(vSoftSpline)
456 >    call deleteSpline(vSoftpSpline)
457  
458 < subroutine useGeometricMixing()
459 <  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
460 <  
461 <  call module_useGeometricMixing()
462 <  return
463 < end subroutine useGeometricMixing
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 >
470 >       call lookupUniformSpline1d(vLJSpline, r, myPot, myDeriv)
471 >      
472 >    else
473 >       ri = 1.0_DP / r
474 >       ri2 = ri*ri
475 >       ri6 = ri2*ri2*ri2
476 >       ri7 = ri6*ri
477 >       ri12 = ri6*ri6
478 >       ri13 = ri12*ri
479 >      
480 >       myPot = 4.0_DP * (ri12 - ri6)
481 >       myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
482 >    endif
483 >
484 >    return
485 >  end subroutine getLJfunc
486 >
487 >  subroutine getSoftFunc(r, myPot, myDeriv)
488 >    
489 >    real(kind=dp), intent(in) :: r
490 >    real(kind=dp), intent(inout) :: myPot, myDeriv
491 >    real(kind=dp) :: ri, ri2, ri6, ri7
492 >    real(kind=dp) :: a, b, c, d, dx
493 >    integer :: j
494 >    
495 >    if (useSplines) then
496 >
497 >       call lookupUniformSpline1d(vSoftSpline, r, myPot, myDeriv)
498 >      
499 >    else
500 >       ri = 1.0_DP / r    
501 >       ri2 = ri*ri
502 >       ri6 = ri2*ri2*ri2
503 >       ri7 = ri6*ri
504 >       myPot = 4.0_DP * (ri6)
505 >       myDeriv = - 24.0_DP * ri7
506 >    endif
507 >
508 >    return
509 >  end subroutine getSoftFunc
510 >
511 >  subroutine setupSplines(rmin, rmax, np)
512 >    real( kind = dp ), intent(in) :: rmin, rmax
513 >    integer, intent(in) :: np
514 >    real( kind = dp ) :: rvals(np), vLJ(np), vLJp(np), vSoft(np), vSoftp(np)
515 >    real( kind = dp ) :: dr, r, ri, ri2, ri6, ri7, ri12, ri13
516 >    integer :: i
517 >
518 >    dr = (rmax-rmin) / float(np-1)
519 >    
520 >    do i = 1, np
521 >       r = rmin + dble(i-1)*dr
522 >       ri = 1.0_DP / r
523 >       ri2 = ri*ri
524 >       ri6 = ri2*ri2*ri2
525 >       ri7 = ri6*ri
526 >       ri12 = ri6*ri6
527 >       ri13 = ri12*ri
528 >
529 >       rvals(i) = r
530 >       vLJ(i) = 4.0_DP * (ri12 - ri6)
531 >       vLJp(i) = 24.0_DP * (ri7 - 2.0_DP * ri13)
532 >
533 >       vSoft(i) = 4.0_DP * (ri6)
534 >       vSoftp(i) = - 24.0_DP * ri7      
535 >    enddo
536 >
537 >    call newSpline(vLJspline, rvals, vLJ,  .true.)
538 >    call newSpline(vLJpspline, rvals, vLJp, .true.)
539 >    call newSpline(vSoftSpline, rvals, vSoft,  .true.)
540 >    call newSpline(vSoftpSpline, rvals, vSoftp, .true.)
541 >
542 >    return
543 >  end subroutine setupSplines
544 > end module lj

Diff Legend

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