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 2062 by tim, Fri Feb 25 21:22:00 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.
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.8 2005-02-25 21:22:00 tim Exp $, $Date: 2005-02-25 21:22:00 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
47  
48 +
49   module lj
50    use atype_module
51    use switcheroo
# Line 20 | Line 63 | module lj
63    integer, parameter :: DP = selected_real_kind(15)
64    
65    type, private :: LjType
66 <     integer :: ident
66 >     integer :: c_ident
67 >     integer :: atid
68       real(kind=dp) :: sigma
69       real(kind=dp) :: epsilon
70 +     logical :: soft_pot
71    end type LjType
72    
73    type(LjType), dimension(:), allocatable :: ParameterMap
# Line 36 | Line 81 | module lj
81       real(kind=dp)  :: tp6
82       real(kind=dp)  :: tp12
83       real(kind=dp)  :: delta
84 +     logical :: soft_pot    
85    end type MixParameters
86    
87    type(MixParameters), dimension(:,:), allocatable :: MixingMap
# Line 56 | Line 102 | contains
102    
103   contains
104  
105 <  subroutine newLJtype(ident, sigma, epsilon, status)
106 <    integer,intent(in) :: ident
105 >  subroutine newLJtype(c_ident, sigma, epsilon, soft_pot, status)
106 >    integer,intent(in) :: c_ident
107      real(kind=dp),intent(in) :: sigma
108      real(kind=dp),intent(in) :: epsilon
109 +    integer, intent(in) :: soft_pot
110      integer,intent(out) :: status
111 <    integer :: nAtypes
111 >    integer :: nATypes, myATID
112 >    integer, pointer :: MatchList(:) => null()
113  
114      status = 0
115      
116 <    !! Be simple-minded and assume that we need a ParameterMap that
69 <    !! is the same size as the total number of atom types
116 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
117  
118      if (.not.allocated(ParameterMap)) then
119        
120 +       !call getMatchingElementList(atypes, "is_LennardJones", .true., &
121 +       !     nLJTypes, MatchList)
122         nAtypes = getSize(atypes)
74 –    
123         if (nAtypes == 0) then
124            status = -1
125            return
# Line 83 | Line 131 | contains
131        
132      end if
133  
134 <    if (ident .gt. size(ParameterMap)) then
134 >    if (myATID .gt. size(ParameterMap)) then
135         status = -1
136         return
137      endif
138      
139      ! set the values for ParameterMap for this atom type:
140  
141 <    ParameterMap(ident)%ident = ident
142 <    ParameterMap(ident)%epsilon = epsilon
143 <    ParameterMap(ident)%sigma = sigma
144 <    
141 >    ParameterMap(myATID)%c_ident = c_ident
142 >    ParameterMap(myATID)%atid = myATID
143 >    ParameterMap(myATID)%epsilon = epsilon
144 >    ParameterMap(myATID)%sigma = sigma
145 >    if (soft_pot == 1) then
146 >      ParameterMap(myATID)%soft_pot = .true.
147 >    else
148 >      ParameterMap(myATID)%soft_pot = .false.
149 >    endif
150    end subroutine newLJtype
151  
152    function getSigma(atid) result (s)
# Line 115 | Line 168 | contains
168      real(kind=dp) :: e
169      
170      if (.not.allocated(ParameterMap)) then
171 <       call handleError("dipole-dipole", "no ParameterMap was present before first call of getEpsilon!")
171 >       call handleError("LJ", "no ParameterMap was present before first call of getEpsilon!")
172         return
173      end if
174      
# Line 148 | Line 201 | contains
201    end subroutine useGeometricMixing
202    
203    subroutine createMixingMap(status)
204 <    integer :: nAtypes
204 >    integer :: nATIDs
205      integer :: status
206      integer :: i
207      integer :: j
# Line 158 | Line 211 | contains
211  
212      status = 0
213      
214 <    nAtypes = size(ParameterMap)
214 >    nATIDs = size(ParameterMap)
215      
216 <    if (nAtypes == 0) then
216 >    if (nATIDs == 0) then
217         status = -1
218         return
219      end if
# Line 171 | Line 224 | contains
224      endif
225      
226      if (.not. allocated(MixingMap)) then
227 <       allocate(MixingMap(nAtypes, nAtypes))
227 >       allocate(MixingMap(nATIDs, nATIDs))
228      endif
229      
230      rcut6 = LJ_rcut**6
231      
232      ! This loops through all atypes, even those that don't support LJ forces.
233 <    do i = 1, nAtypes
233 >    do i = 1, nATIDs
234        
235         Epsilon_i = ParameterMap(i)%epsilon
236         Sigma_i = ParameterMap(i)%sigma
# Line 190 | Line 243 | contains
243         MixingMap(i,i)%epsilon = Epsilon_i          
244         MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
245              (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
246 +       MixingMap(i,i)%soft_pot = ParameterMap(i)%soft_pot
247        
248 <       do j = i + 1, nAtypes
248 >       do j = i + 1, nATIDs
249            
250            Epsilon_j = ParameterMap(j)%epsilon
251            Sigma_j = ParameterMap(j)%sigma
# Line 214 | Line 268 | contains
268            
269            MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
270                 (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
271 +
272 +          MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
273 +
274            
275            MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
276            MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
# Line 221 | Line 278 | contains
278            MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
279            MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
280            MixingMap(j,i)%delta   = MixingMap(i,j)%delta
281 +          MixingMap(j,i)%soft_pot   = MixingMap(i,j)%soft_pot
282            
283         end do
284      end do
285      
286 +    haveMixingMap = .true.
287 +
288    end subroutine createMixingMap
289          
290    subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
# Line 248 | Line 308 | contains
308      real( kind = dp ) :: t6
309      real( kind = dp ) :: t12
310      real( kind = dp ) :: delta
311 +    logical :: soft_pot
312      integer :: id1, id2, localError
313  
314      if (.not.haveMixingMap) then
# Line 264 | Line 325 | contains
325      sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
326      epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
327      delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
328 +    soft_pot =  MixingMap(atid_Row(atom1),atid_Col(atom2))%soft_pot
329   #else
330      sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
331      epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
332      delta    = MixingMap(atid(atom1),atid(atom2))%delta
333 +    soft_pot    = MixingMap(atid(atom1),atid(atom2))%soft_pot
334   #endif
335  
336      r6 = r2 * r2 * r2
337      
338      t6  = sigma6/ r6
339      t12 = t6 * t6    
277 –  
278 –    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
279 –    if (LJ_do_shift) then
280 –       pot_temp = pot_temp + delta
281 –    endif
340  
341 <    vpair = vpair + pot_temp
342 <      
343 <    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
341 >    if (soft_pot) then
342 >      pot_temp = 4.0E0_DP * epsilon * t12
343 >      if (LJ_do_shift) then
344 >         pot_temp = pot_temp + delta
345 >      endif
346 >  
347 >      vpair = vpair + pot_temp
348 >        
349 >      dudr = sw * 24.0E0_DP * epsilon * (-2.0E0_DP)*t12 / rij
350 >
351 >    else
352 >      pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
353 >      if (LJ_do_shift) then
354 >         pot_temp = pot_temp + delta
355 >      endif
356 >  
357 >      vpair = vpair + pot_temp
358 >        
359 >      dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
360 >    endif
361        
362      drdx = d(1) / rij
363      drdy = d(2) / rij
# Line 343 | Line 418 | end module lj
418    !! Calculates the mixing for sigma or epslon
419      
420   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 –  call module_newLJtype(ident, sigma, epsilon, status)
356 –  
357 – end subroutine newLJtype
358 –
359 – subroutine useGeometricMixing()
360 –  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
361 –  
362 –  call module_useGeometricMixing()
363 –  return
364 – end subroutine useGeometricMixing

Diff Legend

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