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 1624 by chuckv, Thu Oct 21 15:25:30 2004 UTC vs.
Revision 2461 by gezelter, Mon Nov 21 22:59:02 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.2 2004-10-21 15:25:30 chuckv Exp $, $Date: 2004-10-21 15:25:30 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
46 > !! @version $Id: LJ.F90,v 1.19 2005-11-21 22:59:01 gezelter Exp $, $Date: 2005-11-21 22:59:01 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
47  
48 +
49   module lj
50    use atype_module
9  use switcheroo
51    use vector_class
52    use simulation
53 +  use status
54   #ifdef IS_MPI
55    use mpiSimulation
56   #endif
# Line 16 | Line 58 | module lj
58  
59    implicit none
60    PRIVATE
61 <  
61 > #define __FORTRAN90
62 > #include "UseTheForce/DarkSide/fInteractionMap.h"
63 >
64    integer, parameter :: DP = selected_real_kind(15)
65  
66 < #define __FORTRAN90
67 < #include "UseTheForce/fForceField.h"
66 >  logical, save :: useGeometricDistanceMixing = .false.
67 >  logical, save :: haveMixingMap = .false.
68  
69 <  integer, save :: LJ_Mixing_Policy
70 <  real(kind=DP), save :: LJ_rcut
71 <  logical, save :: havePolicy = .false.
72 <  logical, save :: haveCut = .false.
73 <  logical, save :: LJ_do_shift = .false.
74 <  
75 <  !! Logical has lj force field module been initialized?
76 <  
77 <  logical, save :: LJ_FF_initialized = .false.
78 <  
79 <  !! Public methods and data
80 <  public :: init_LJ_FF
81 <  public :: setCutoffLJ
82 <  public :: do_lj_pair
69 >  real(kind=DP), save :: defaultCutoff = 0.0_DP
70 >  logical, save :: defaultShift = .false.
71 >  logical, save :: haveDefaultCutoff = .false.
72 >
73 >
74 >  type, private :: LJtype
75 >     integer       :: atid
76 >     real(kind=dp) :: sigma
77 >     real(kind=dp) :: epsilon
78 >     logical       :: isSoftCore = .false.
79 >  end type LJtype
80 >
81 >  type, private :: LJList
82 >     integer               :: Nljtypes = 0
83 >     integer               :: currentLJtype = 0
84 >     type(LJtype), pointer :: LJtypes(:)      => null()
85 >     integer, pointer      :: atidToLJtype(:) => null()
86 >  end type LJList
87 >
88 >  type(LJList), save :: LJMap
89 >
90 >  type :: MixParameters
91 >     real(kind=DP) :: sigma
92 >     real(kind=DP) :: epsilon
93 >     real(kind=dp) :: sigma6
94 >     real(kind=dp) :: rCut
95 >     real(kind=dp) :: delta
96 >     logical       :: rCutWasSet = .false.
97 >     logical       :: shiftedPot
98 >     logical       :: isSoftCore = .false.
99 >  end type MixParameters
100 >
101 >  type(MixParameters), dimension(:,:), allocatable :: MixingMap
102 >
103    public :: newLJtype
104 <  
105 <  !! structure for lj type parameters
106 <  type, private :: ljType
107 <    integer :: lj_ident
108 <    real(kind=dp) :: lj_sigma
109 <    real(kind=dp) :: lj_epsilon
110 <  end type ljType
47 <  
48 <  !! List of lj type parameters
49 <  type, private :: ljTypeList
50 <    integer  :: n_lj_types = 0
51 <    integer  :: currentAddition = 0
52 <    type(ljType), pointer :: ljParams(:) => null()
53 <  end type ljTypeList
54 <  
55 <  !! The list of lj Parameters
56 <  type (ljTypeList), save :: ljParameterList
57 <  
58 <  
59 <  type :: lj_mixed_params
60 <     !! Lennard-Jones epsilon
61 <     real ( kind = dp )  :: epsilon = 0.0_dp
62 <     !! Lennard-Jones Sigma
63 <     real ( kind = dp )  :: sigma = 0.0_dp
64 <     !! Lennard-Jones Sigma to sixth
65 <     real ( kind = dp )  :: sigma6 = 0.0_dp
66 <     !!
67 <     real ( kind = dp )  :: tp6
68 <     real ( kind = dp )  :: tp12
69 <     real ( kind = dp )  :: delta  = 0.0_dp
70 <  end type lj_mixed_params
71 <  
72 <  type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
73 <  
74 <  
75 <  
104 >  public :: setLJDefaultCutoff
105 >  public :: getSigma
106 >  public :: getEpsilon
107 >  public :: useGeometricMixing
108 >  public :: do_lj_pair
109 >  public :: destroyLJtypes
110 >
111   contains
112  
113 <  subroutine newLJtype(ident,lj_sigma,lj_epsilon,status)
114 <    integer,intent(in) :: ident
115 <    real(kind=dp),intent(in) :: lj_sigma
116 <    real(kind=dp),intent(in) :: lj_epsilon
113 >  subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
114 >    integer,intent(in) :: c_ident
115 >    real(kind=dp),intent(in) :: sigma
116 >    real(kind=dp),intent(in) :: epsilon
117 >    integer, intent(in) :: isSoftCore
118      integer,intent(out) :: status
119 <    
120 <    integer,pointer                        :: Matchlist(:) => null()
119 >    integer :: nLJTypes, ntypes, myATID
120 >    integer, pointer :: MatchList(:) => null()
121      integer :: current
122 <    integer :: nAtypes
122 >
123      status = 0
124 <    
125 <        !! Assume that atypes has already been set and get the total number of types in atypes
90 <  
91 <  
124 >    ! check to see if this is the first time into this routine...
125 >    if (.not.associated(LJMap%LJtypes)) then
126  
127 <    ! check to see if this is the first time into
128 <    if (.not.associated(ljParameterList%ljParams)) then
129 <       call getMatchingElementList(atypes, "is_lj", .true., nAtypes, MatchList)
130 <       ljParameterList%n_lj_types = nAtypes
131 <       if (nAtypes == 0) then
132 <         status = -1
133 <         return
134 <       end if
135 <       allocate(ljParameterList%ljParams(nAtypes))
127 >       call getMatchingElementList(atypes, "is_LennardJones", .true., &
128 >            nLJTypes, MatchList)
129 >      
130 >       LJMap%nLJtypes =  nLJTypes
131 >
132 >       allocate(LJMap%LJtypes(nLJTypes))
133 >
134 >       ntypes = getSize(atypes)
135 >
136 >       allocate(LJMap%atidToLJtype(ntypes))
137      end if
138  
139 <    ljParameterList%currentAddition = ljParameterList%currentAddition + 1
140 <    current = ljParameterList%currentAddition
141 <    
142 <    ! set the values for ljParameterList
143 <    ljParameterList%ljParams(current)%lj_ident = ident
144 <    ljParameterList%ljParams(current)%lj_epsilon = lj_epsilon
145 <    ljParameterList%ljParams(current)%lj_sigma = lj_sigma
146 <    
147 <  end subroutine newLJtype
148 <  
114 <  subroutine init_LJ_FF(mix_Policy, status)
115 <    integer, intent(in) :: mix_Policy
116 <    integer, intent(out) :: status
117 <    integer :: myStatus
118 <    
119 <    if (mix_Policy == LB_MIXING_RULE) then
120 <       LJ_Mixing_Policy = LB_MIXING_RULE
139 >    LJMap%currentLJtype = LJMap%currentLJtype + 1
140 >    current = LJMap%currentLJtype
141 >
142 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
143 >    LJMap%atidToLJtype(myATID)        = current
144 >    LJMap%LJtypes(current)%atid       = myATID
145 >    LJMap%LJtypes(current)%sigma      = sigma
146 >    LJMap%LJtypes(current)%epsilon    = epsilon
147 >    if (isSoftCore .eq. 1) then
148 >       LJMap%LJtypes(current)%isSoftCore = .true.
149      else
150 <       if (mix_Policy == EXPLICIT_MIXING_RULE) then
123 <          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
124 <       else
125 <          write(*,*) 'Unknown Mixing Policy!'
126 <          status = -1
127 <          return
128 <       endif
150 >       LJMap%LJtypes(current)%isSoftCore = .false.
151      endif
152 +  end subroutine newLJtype
153  
154 <    havePolicy = .true.
154 >  subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
155 >    real(kind=dp), intent(in) :: thisRcut
156 >    logical, intent(in) :: shiftedPot
157 >    defaultCutoff = thisRcut
158 >    defaultShift = shiftedPot
159 >    haveDefaultCutoff = .true.
160 >    !we only want to build LJ Mixing map if LJ is being used.
161 >    if(LJMap%nLJTypes /= 0) then
162 >       call createMixingMap()
163 >    end if
164 >  end subroutine setLJDefaultCutoff
165  
166 <    if (haveCut) then
167 <       status = 0
168 <       call createMixingList(myStatus)
169 <       if (myStatus /= 0) then
170 <          status = -1
171 <          return
172 <       end if
173 <      
141 <       LJ_FF_initialized = .true.
166 >  function getSigma(atid) result (s)
167 >    integer, intent(in) :: atid
168 >    integer :: ljt1
169 >    real(kind=dp) :: s
170 >
171 >    if (LJMap%currentLJtype == 0) then
172 >       call handleError("LJ", "No members in LJMap")
173 >       return
174      end if
143  
144  end subroutine init_LJ_FF
145  
146  subroutine setCutoffLJ(rcut, do_shift, status)
147    logical, intent(in):: do_shift
148    integer :: status, myStatus
149    real(kind=dp) :: rcut
175  
176 < #define __FORTRAN90
177 < #include "UseTheForce/fSwitchingFunction.h"
176 >    ljt1 = LJMap%atidToLJtype(atid)
177 >    s = LJMap%LJtypes(ljt1)%sigma
178  
179 <    status = 0
179 >  end function getSigma
180  
181 <    LJ_rcut = rcut
182 <    LJ_do_shift = do_shift
183 <    call set_switch(LJ_SWITCH, rcut, rcut)
184 <    haveCut = .true.
181 >  function getEpsilon(atid) result (e)
182 >    integer, intent(in) :: atid
183 >    integer :: ljt1
184 >    real(kind=dp) :: e
185  
186 <    if (havePolicy) then
187 <       status = 0
188 <       call createMixingList(myStatus)
189 <       if (myStatus /= 0) then
190 <          status = -1
191 <          return
192 <       end if
193 <      
194 <       LJ_FF_initialized = .true.
195 <    end if    
196 <    
186 >    if (LJMap%currentLJtype == 0) then
187 >       call handleError("LJ", "No members in LJMap")
188 >       return
189 >    end if
190 >
191 >    ljt1 = LJMap%atidToLJtype(atid)
192 >    e = LJMap%LJtypes(ljt1)%epsilon
193 >
194 >  end function getEpsilon
195 >
196 >  subroutine useGeometricMixing()
197 >    useGeometricDistanceMixing = .true.
198 >    haveMixingMap = .false.
199      return
200 <  end subroutine setCutoffLJ
201 <  
202 <  subroutine createMixingList(status)
203 <    integer :: nAtypes
204 <    integer :: status
205 <    integer :: i
206 <    integer :: j
207 <    real ( kind = dp ) :: mySigma_i,mySigma_j
208 <    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
209 <    real ( kind = dp ) :: rcut6
183 <    logical :: I_isLJ, J_isLJ
184 <    status = 0
185 <    
186 <    ! we only allocate this array to the number of lj_atypes
187 <    nAtypes = size(ljParameterList%ljParams)
188 <    if (nAtypes == 0) then
189 <       status = -1
200 >  end subroutine useGeometricMixing
201 >
202 >  subroutine createMixingMap()
203 >    integer :: nLJtypes, i, j
204 >    real ( kind = dp ) :: s1, s2, e1, e2
205 >    real ( kind = dp ) :: rcut6, tp6, tp12
206 >    logical :: isSoftCore1, isSoftCore2, doShift
207 >
208 >    if (LJMap%currentLJtype == 0) then
209 >       call handleError("LJ", "No members in LJMap")
210         return
211      end if
212 <        
213 <    if (.not. associated(ljMixed)) then
214 <       allocate(ljMixed(nAtypes, nAtypes))
212 >
213 >    nLJtypes = LJMap%nLJtypes
214 >
215 >    if (.not. allocated(MixingMap)) then
216 >       allocate(MixingMap(nLJtypes, nLJtypes))
217      endif
218  
219 <    rcut6 = LJ_rcut**6
219 >    do i = 1, nLJtypes
220  
221 < ! This loops through all atypes, even those that don't support LJ forces.
222 <    do i = 1, nAtypes
221 >       s1 = LJMap%LJtypes(i)%sigma
222 >       e1 = LJMap%LJtypes(i)%epsilon
223 >       isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
224  
225 <          myEpsilon_i = ljParameterList%ljParams(i)%lj_epsilon
203 <          mySigma_i = ljParameterList%ljParams(i)%lj_sigma
225 >       do j = i, nLJtypes
226            
227 <          ! do self mixing rule
228 <          ljMixed(i,i)%sigma   = mySigma_i
227 >          s2 = LJMap%LJtypes(j)%sigma
228 >          e2 = LJMap%LJtypes(j)%epsilon
229 >          isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
230            
231 <          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
231 >          MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
232 >
233 >          ! only the distance parameter uses different mixing policies
234 >          if (useGeometricDistanceMixing) then
235 >             MixingMap(i,j)%sigma = dsqrt(s1 * s2)
236 >          else
237 >             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
238 >          endif
239            
240 <          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
211 <          
212 <          ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
213 <          
214 <          
215 <          ljMixed(i,i)%epsilon = myEpsilon_i
216 <          
217 <          ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
218 <            (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
219 <          
220 <          do j = i + 1, nAtypes
240 >          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
241  
242 <                myEpsilon_j = ljParameterList%ljParams(j)%lj_epsilon
223 <                mySigma_j = ljParameterList%ljParams(j)%lj_sigma
242 >          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
243  
244 <                          
245 <                ljMixed(i,j)%sigma  =  &
246 <                     calcLJMix("sigma",mySigma_i, &
247 <                     mySigma_j)
248 <                
249 <                ljMixed(i,j)%sigma6 = &
250 <                     (ljMixed(i,j)%sigma)**6
251 <                
252 <                
253 <                ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
254 <                
255 <                ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
256 <                
257 <                
258 <                ljMixed(i,j)%epsilon = &
259 <                     calcLJMix("epsilon",myEpsilon_i, &
260 <                     myEpsilon_j)
261 <                
262 <                ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
263 <                     (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
264 <                
265 <                
266 <                ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
267 <                ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
249 <                ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
250 <                ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
251 <                ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
252 <                ljMixed(j,i)%delta   = ljMixed(i,j)%delta
253 <          
254 <          end do
255 <    end do
244 >          if (haveDefaultCutoff) then
245 >             rcut6 = defaultCutoff**6
246 >             tp6    = MixingMap(i,j)%sigma6/rcut6
247 >             tp12    = tp6**2          
248 >             MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
249 >             MixingMap(i,j)%shiftedPot = defaultShift
250 >          else
251 >             MixingMap(i,j)%delta = 0.0_DP
252 >             MixingMap(i,j)%shiftedPot = defaultShift
253 >          endif          
254 >
255 >          if (i.ne.j) then
256 >             MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
257 >             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
258 >             MixingMap(j,i)%sigma6     = MixingMap(i,j)%sigma6
259 >             MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
260 >             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
261 >             MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
262 >             MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
263 >             MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
264 >          endif
265 >
266 >       enddo
267 >    enddo
268      
269 <  end subroutine createMixingList
269 >    haveMixingMap = .true.
270 >    
271 >  end subroutine createMixingMap
272    
273 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
273 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
274         pot, f, do_pot)
275 <
275 >    
276      integer, intent(in) ::  atom1, atom2
277 <    real( kind = dp ), intent(in) :: rij, r2
277 >    integer :: atid1, atid2, ljt1, ljt2
278 >    real( kind = dp ), intent(in) :: rij, r2, rcut
279      real( kind = dp ) :: pot, sw, vpair
280      real( kind = dp ), dimension(3,nLocal) :: f    
281      real( kind = dp ), intent(in), dimension(3) :: d
# Line 273 | Line 288 | contains
288      real( kind = dp ) :: pot_temp, dudr
289      real( kind = dp ) :: sigma6
290      real( kind = dp ) :: epsilon
291 <    real( kind = dp ) :: r6
292 <    real( kind = dp ) :: t6
293 <    real( kind = dp ) :: t12
291 >    real( kind = dp ) :: r6, rc6
292 >    real( kind = dp ) :: t6, tp6
293 >    real( kind = dp ) :: t12, tp12
294      real( kind = dp ) :: delta
295 <    integer :: id1, id2
295 >    logical :: isSoftCore, shiftedPot
296 >    integer :: id1, id2, localError
297  
298 +    if (.not.haveMixingMap) then
299 +       call createMixingMap()
300 +    endif
301 +
302      ! Look up the correct parameters in the mixing matrix
303   #ifdef IS_MPI
304 <    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
305 <    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
286 <    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
304 >    atid1 = atid_Row(atom1)
305 >    atid2 = atid_Col(atom2)
306   #else
307 <    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
308 <    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
290 <    delta    = ljMixed(atid(atom1),atid(atom2))%delta
307 >    atid1 = atid(atom1)
308 >    atid2 = atid(atom2)
309   #endif
310  
311 +    ljt1 = LJMap%atidToLJtype(atid1)
312 +    ljt2 = LJMap%atidToLJtype(atid2)
313 +
314 +    sigma6     = MixingMap(ljt1,ljt2)%sigma6
315 +    epsilon    = MixingMap(ljt1,ljt2)%epsilon
316 +    delta      = MixingMap(ljt1,ljt2)%delta
317 +    isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
318 +    shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
319 +
320      r6 = r2 * r2 * r2
321 <    
321 >
322      t6  = sigma6/ r6
323      t12 = t6 * t6    
297  
298    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
299    if (LJ_do_shift) then
300       pot_temp = pot_temp + delta
301    endif
324  
325 <    vpair = vpair + pot_temp
325 >    if (isSoftCore) then
326        
327 <    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
327 >       pot_temp = 4.0E0_DP * epsilon * t6
328 >       if (shiftedPot) then
329 >          rc6 = rcut**6
330 >          tp6 = sigma6 / rc6
331 >          delta =-4.0_DP*epsilon*(tp6)
332 >          pot_temp = pot_temp + delta
333 >       endif
334        
335 +       vpair = vpair + pot_temp
336 +      
337 +       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
338 +      
339 +    else
340 +       pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
341 +       if (shiftedPot) then
342 +          rc6 = rcut**6
343 +          tp6 = sigma6 / rc6
344 +          tp12 = tp6*tp6
345 +          delta =-4.0_DP*epsilon*(tp12 - tp6)
346 +          pot_temp = pot_temp + delta
347 +       endif
348 +      
349 +       vpair = vpair + pot_temp
350 +      
351 +       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
352 +    endif
353 +
354      drdx = d(1) / rij
355      drdy = d(2) / rij
356      drdz = d(3) / rij
357 <      
357 >
358      fx = dudr * drdx
359      fy = dudr * drdy
360      fz = dudr * drdz
361 <    
315 <      
361 >
362   #ifdef IS_MPI
363      if (do_pot) then
364 <       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
365 <       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
364 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
365 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
366      endif
367 <    
367 >
368      f_Row(1,atom1) = f_Row(1,atom1) + fx
369      f_Row(2,atom1) = f_Row(2,atom1) + fy
370      f_Row(3,atom1) = f_Row(3,atom1) + fz
371 <    
371 >
372      f_Col(1,atom2) = f_Col(1,atom2) - fx
373      f_Col(2,atom2) = f_Col(2,atom2) - fy
374      f_Col(3,atom2) = f_Col(3,atom2) - fz      
375 <    
375 >
376   #else
377      if (do_pot) pot = pot + sw*pot_temp
378  
379      f(1,atom1) = f(1,atom1) + fx
380      f(2,atom1) = f(2,atom1) + fy
381      f(3,atom1) = f(3,atom1) + fz
382 <    
382 >
383      f(1,atom2) = f(1,atom2) - fx
384      f(2,atom2) = f(2,atom2) - fy
385      f(3,atom2) = f(3,atom2) - fz
386   #endif
387 <        
387 >
388   #ifdef IS_MPI
389      id1 = AtomRowToGlobal(atom1)
390      id2 = AtomColToGlobal(atom2)
# Line 348 | Line 394 | contains
394   #endif
395  
396      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
397 <      
397 >
398         fpair(1) = fpair(1) + fx
399         fpair(2) = fpair(2) + fy
400         fpair(3) = fpair(3) + fz
# Line 356 | Line 402 | contains
402      endif
403  
404      return    
405 <    
405 >
406    end subroutine do_lj_pair
361  
362  
363  !! Calculates the mixing for sigma or epslon
364  
365  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
366    character(len=*) :: thisParam
367    real(kind = dp)  :: param1
368    real(kind = dp)  :: param2
369    real(kind = dp ) :: myMixParam
407  
408 <    integer, optional :: status  
408 >  subroutine destroyLJTypes()
409  
410 <    myMixParam = 0.0_dp
410 >    LJMap%nLJtypes = 0
411 >    LJMap%currentLJtype = 0
412      
413 <    if (present(status)) status = 0
414 <    select case (LJ_Mixing_Policy)
415 <    case (1)
416 <       select case (thisParam)
417 <       case ("sigma")
418 <          myMixParam = 0.5_dp * (param1 + param2)
419 <       case ("epsilon")
420 <          myMixParam = sqrt(param1 * param2)
421 <       case default
422 <          status = -1
423 <       end select
424 <    case default
387 <       status = -1
388 <    end select
389 <  end function calcLJMix
390 <  
391 < end module lj
413 >    if (associated(LJMap%LJtypes)) then
414 >       deallocate(LJMap%LJtypes)
415 >       LJMap%LJtypes => null()
416 >    end if
417 >    
418 >    if (associated(LJMap%atidToLJtype)) then
419 >       deallocate(LJMap%atidToLJtype)
420 >       LJMap%atidToLJtype => null()
421 >    end if
422 >    
423 >    haveMixingMap = .false.
424 >  end subroutine destroyLJTypes
425  
426 < subroutine newLJtype(ident,lj_sigma,lj_epsilon,status)
394 <    use lj, ONLY : module_newLJtype => newLJtype
395 <    integer, parameter :: DP = selected_real_kind(15)
396 <    integer,intent(inout) :: ident
397 <    real(kind=dp),intent(inout) :: lj_sigma
398 <    real(kind=dp),intent(inout) :: lj_epsilon
399 <    integer,intent(inout) :: status
400 <
401 <    call module_newLJtype(ident,lj_sigma,lj_epsilon,status)
402 <
403 < end subroutine newLJtype
404 <
426 > end module lj

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines