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 2062 by tim, Fri Feb 25 21:22:00 2005 UTC vs.
Revision 2189 by chuckv, Wed Apr 13 20:36:45 2005 UTC

# Line 43 | Line 43
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.8 2005-02-25 21:22:00 tim Exp $, $Date: 2005-02-25 21:22:00 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
46 > !! @version $Id: LJ.F90,v 1.10 2005-04-13 20:36:45 chuckv Exp $, $Date: 2005-04-13 20:36:45 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $
47  
48  
49   module lj
# Line 99 | Line 99 | module lj
99    public :: newLJtype  
100    public :: getSigma
101    public :: getEpsilon
102 +  public :: destroyLJTypes
103    
104   contains
105  
# Line 124 | Line 125 | contains
125            status = -1
126            return
127         end if
128 <      
128 >        
129         if (.not. allocated(ParameterMap)) then
130            allocate(ParameterMap(nAtypes))
131         endif
132 <      
132 >      
133      end if
134  
135      if (myATID .gt. size(ParameterMap)) then
# Line 208 | Line 209 | contains
209      real ( kind = dp ) :: Sigma_i, Sigma_j
210      real ( kind = dp ) :: Epsilon_i, Epsilon_j
211      real ( kind = dp ) :: rcut6
212 +    logical :: i_is_LJ, j_is_LJ
213  
214      status = 0
215 +
216 +    if (.not. allocated(ParameterMap)) then
217 +       call handleError("LJ", "no ParameterMap was present before call of createMixingMap!")
218 +       status = -1
219 +       return
220 +    endif
221      
222      nATIDs = size(ParameterMap)
223      
# Line 218 | Line 226 | contains
226         return
227      end if
228  
229 +    if (.not. allocated(MixingMap)) then
230 +       allocate(MixingMap(nATIDs, nATIDs))
231 +    endif
232 +
233      if (.not.have_rcut) then
234         status = -1
235         return
236      endif
237 <    
226 <    if (.not. allocated(MixingMap)) then
227 <       allocate(MixingMap(nATIDs, nATIDs))
228 <    endif
229 <    
237 >        
238      rcut6 = LJ_rcut**6
239      
240      ! This loops through all atypes, even those that don't support LJ forces.
241      do i = 1, nATIDs
242 <      
243 <       Epsilon_i = ParameterMap(i)%epsilon
244 <       Sigma_i = ParameterMap(i)%sigma
245 <      
238 <       ! do self mixing rule
239 <       MixingMap(i,i)%sigma   = Sigma_i          
240 <       MixingMap(i,i)%sigma6  = Sigma_i ** 6          
241 <       MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
242 <       MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
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, nATIDs
242 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
243 >       if (i_is_LJ) then
244 >          Epsilon_i = ParameterMap(i)%epsilon
245 >          Sigma_i = ParameterMap(i)%sigma
246            
247 <          Epsilon_j = ParameterMap(j)%epsilon
248 <          Sigma_j = ParameterMap(j)%sigma
247 >          ! do self mixing rule
248 >          MixingMap(i,i)%sigma   = Sigma_i          
249 >          MixingMap(i,i)%sigma6  = Sigma_i ** 6          
250 >          MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
251 >          MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
252 >          MixingMap(i,i)%epsilon = Epsilon_i          
253 >          MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
254 >               (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
255 >          MixingMap(i,i)%soft_pot = ParameterMap(i)%soft_pot
256            
257 <          ! only the distance parameter uses different mixing policies
258 <          if (useGeometricDistanceMixing) then
255 <             ! only for OPLS as far as we can tell
256 <             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
257 <          else
258 <             ! everyone else
259 <             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
260 <          endif
257 >          do j = i + 1, nATIDs
258 >             call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
259            
260 <          ! energy parameter is always geometric mean:
261 <          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
262 <                    
263 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
264 <          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
265 <          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
266 <          
267 <          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
268 <               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
269 <
270 <          MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
271 <
272 <          
273 <          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
274 <          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
275 <          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
276 <          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
277 <          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
278 <          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
279 <          MixingMap(j,i)%soft_pot   = MixingMap(i,j)%soft_pot
280 <          
281 <       end do
260 >             if (j_is_LJ) then
261 >                Epsilon_j = ParameterMap(j)%epsilon
262 >                Sigma_j = ParameterMap(j)%sigma
263 >                
264 >                ! only the distance parameter uses different mixing policies
265 >                if (useGeometricDistanceMixing) then
266 >                   ! only for OPLS as far as we can tell
267 >                   MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
268 >                else
269 >                   ! everyone else
270 >                   MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
271 >                endif
272 >                
273 >                ! energy parameter is always geometric mean:
274 >                MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
275 >                
276 >                MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
277 >                MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
278 >                MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
279 >                
280 >                MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
281 >                     (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
282 >                
283 >                MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
284 >                
285 >                
286 >                MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
287 >                MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
288 >                MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
289 >                MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
290 >                MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
291 >                MixingMap(j,i)%delta   = MixingMap(i,j)%delta
292 >                MixingMap(j,i)%soft_pot   = MixingMap(i,j)%soft_pot
293 >             endif
294 >          end do
295 >       endif
296      end do
297      
298      haveMixingMap = .true.
# Line 414 | Line 426 | contains
426      
427    end subroutine do_lj_pair
428    
429 <  
430 <  !! Calculates the mixing for sigma or epslon
429 >  subroutine destroyLJTypes()
430 >    if(allocated(ParameterMap)) deallocate(ParameterMap)
431 >    if(allocated(MixingMap)) deallocate(MixingMap)
432 >    haveMixingMap = .false.
433 >  end subroutine destroyLJTypes
434 >
435      
436   end module lj

Diff Legend

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