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 2271 by gezelter, Tue Aug 9 22:33:48 2005 UTC vs.
Revision 2497 by chuckv, Wed Dec 7 19:58:18 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.13 2005-08-09 22:33:48 gezelter Exp $, $Date: 2005-08-09 22:33:48 $, $Name: not supported by cvs2svn $, $Revision: 1.13 $
46 > !! @version $Id: LJ.F90,v 1.20 2005-12-07 19:58:18 chuckv Exp $, $Date: 2005-12-07 19:58:18 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
47  
48  
49   module lj
# Line 51 | Line 51 | module lj
51    use vector_class
52    use simulation
53    use status
54 +  use fForceOptions
55   #ifdef IS_MPI
56    use mpiSimulation
57   #endif
# Line 58 | Line 59 | module lj
59  
60    implicit none
61    PRIVATE
62 + #define __FORTRAN90
63 + #include "UseTheForce/DarkSide/fInteractionMap.h"
64  
65    integer, parameter :: DP = selected_real_kind(15)
66  
# Line 77 | Line 80 | module lj
80    end type LJtype
81  
82    type, private :: LJList
83 <     integer               :: nLJtypes = 0
83 >     integer               :: Nljtypes = 0
84       integer               :: currentLJtype = 0
85       type(LJtype), pointer :: LJtypes(:)      => null()
86       integer, pointer      :: atidToLJtype(:) => null()
# Line 100 | Line 103 | module lj
103  
104    public :: newLJtype
105    public :: setLJDefaultCutoff
103 –  public :: setLJUniformCutoff
104 –  public :: setLJCutoffByTypes
106    public :: getSigma
107    public :: getEpsilon
107 –  public :: useGeometricMixing
108    public :: do_lj_pair
109    public :: destroyLJtypes
110  
# Line 157 | Line 157 | contains
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  
162 –  subroutine setLJUniformCutoff(thisRcut, shiftedPot)
163 –    real(kind=dp), intent(in) :: thisRcut
164 –    logical, intent(in) :: shiftedPot
165 –    integer :: nLJtypes, i, j
166 –
167 –    if (LJMap%currentLJtype == 0) then
168 –       call handleError("LJ", "No members in LJMap")
169 –       return
170 –    end if
171 –
172 –    nLJtypes = LJMap%nLJtypes
173 –    if (.not. allocated(MixingMap)) then
174 –       allocate(MixingMap(nLJtypes, nLJtypes))
175 –    endif
176 –
177 –    do i = 1, nLJtypes
178 –       do j = 1, nLJtypes
179 –          MixingMap(i,j)%rCut = thisRcut
180 –          MixingMap(i,j)%shiftedPot = shiftedPot
181 –          MixingMap(i,j)%rCutWasSet = .true.
182 –       enddo
183 –    enddo
184 –    call createMixingMap()
185 –  end subroutine setLJUniformCutoff
186 –
187 –  subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
188 –    integer, intent(in) :: atid1, atid2
189 –    real(kind=dp), intent(in) :: thisRcut
190 –    logical, intent(in) :: shiftedPot
191 –    integer :: nLJtypes, ljt1, ljt2
192 –
193 –    if (LJMap%currentLJtype == 0) then
194 –       call handleError("LJ", "No members in LJMap")
195 –       return
196 –    end if
197 –
198 –    nLJtypes = LJMap%nLJtypes
199 –    if (.not. allocated(MixingMap)) then
200 –       allocate(MixingMap(nLJtypes, nLJtypes))
201 –    endif
202 –
203 –    ljt1 = LJMap%atidToLJtype(atid1)
204 –    ljt2 = LJMap%atidToLJtype(atid2)
205 –
206 –    MixingMap(ljt1,ljt2)%rCut = thisRcut
207 –    MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
208 –    MixingMap(ljt1,ljt2)%rCutWasSet = .true.
209 –    MixingMap(ljt2,ljt1)%rCut = thisRcut
210 –    MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
211 –    MixingMap(ljt2,ljt1)%rCutWasSet = .true.
212 –
213 –    call createMixingMap()
214 –  end subroutine setLJCutoffByTypes
215 –
166    function getSigma(atid) result (s)
167      integer, intent(in) :: atid
168      integer :: ljt1
# Line 243 | Line 193 | contains
193  
194    end function getEpsilon
195  
246 –  subroutine useGeometricMixing()
247 –    useGeometricDistanceMixing = .true.
248 –    haveMixingMap = .false.
249 –    return
250 –  end subroutine useGeometricMixing
251 –
196    subroutine createMixingMap()
197      integer :: nLJtypes, i, j
198      real ( kind = dp ) :: s1, s2, e1, e2
199      real ( kind = dp ) :: rcut6, tp6, tp12
200 <    logical :: isSoftCore1, isSoftCore2
200 >    logical :: isSoftCore1, isSoftCore2, doShift
201  
202      if (LJMap%currentLJtype == 0) then
203         call handleError("LJ", "No members in LJMap")
# Line 266 | Line 210 | contains
210         allocate(MixingMap(nLJtypes, nLJtypes))
211      endif
212  
213 +    useGeometricDistanceMixing = usesGeometricDistanceMixing()
214      do i = 1, nLJtypes
215  
216         s1 = LJMap%LJtypes(i)%sigma
# Line 291 | Line 236 | contains
236  
237            MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
238  
239 <          if (MixingMap(i,j)%rCutWasSet) then
240 <             rcut6 = (MixingMap(i,j)%rcut)**6
239 >          if (haveDefaultCutoff) then
240 >             rcut6 = defaultCutoff**6
241 >             tp6    = MixingMap(i,j)%sigma6/rcut6
242 >             tp12    = tp6**2          
243 >             MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
244 >             MixingMap(i,j)%shiftedPot = defaultShift
245            else
246 <             if (haveDefaultCutoff) then
247 <                rcut6 = defaultCutoff**6
248 <             else
249 <                call handleError("LJ", "No specified or default cutoff value!")
250 <             endif
246 >             MixingMap(i,j)%delta = 0.0_DP
247 >             MixingMap(i,j)%shiftedPot = defaultShift
248 >          endif          
249 >
250 >          if (i.ne.j) then
251 >             MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
252 >             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
253 >             MixingMap(j,i)%sigma6     = MixingMap(i,j)%sigma6
254 >             MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
255 >             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
256 >             MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
257 >             MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
258 >             MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
259            endif
303 –          
304 –          tp6    = MixingMap(i,j)%sigma6/rcut6
305 –          tp12    = tp6**2          
260  
307 –          MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
261         enddo
262      enddo
263      
# Line 312 | Line 265 | contains
265      
266    end subroutine createMixingMap
267    
268 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
268 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
269         pot, f, do_pot)
270 <
270 >    
271      integer, intent(in) ::  atom1, atom2
272      integer :: atid1, atid2, ljt1, ljt2
273 <    real( kind = dp ), intent(in) :: rij, r2
273 >    real( kind = dp ), intent(in) :: rij, r2, rcut
274      real( kind = dp ) :: pot, sw, vpair
275      real( kind = dp ), dimension(3,nLocal) :: f    
276      real( kind = dp ), intent(in), dimension(3) :: d
# Line 330 | Line 283 | contains
283      real( kind = dp ) :: pot_temp, dudr
284      real( kind = dp ) :: sigma6
285      real( kind = dp ) :: epsilon
286 <    real( kind = dp ) :: r6
287 <    real( kind = dp ) :: t6
288 <    real( kind = dp ) :: t12
286 >    real( kind = dp ) :: r6, rc6
287 >    real( kind = dp ) :: t6, tp6
288 >    real( kind = dp ) :: t12, tp12
289      real( kind = dp ) :: delta
290      logical :: isSoftCore, shiftedPot
291      integer :: id1, id2, localError
# Line 368 | Line 321 | contains
321        
322         pot_temp = 4.0E0_DP * epsilon * t6
323         if (shiftedPot) then
324 +          rc6 = rcut**6
325 +          tp6 = sigma6 / rc6
326 +          delta =-4.0_DP*epsilon*(tp6)
327            pot_temp = pot_temp + delta
328         endif
329        
# Line 378 | Line 334 | contains
334      else
335         pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
336         if (shiftedPot) then
337 +          rc6 = rcut**6
338 +          tp6 = sigma6 / rc6
339 +          tp12 = tp6*tp6
340 +          delta =-4.0_DP*epsilon*(tp12 - tp6)
341            pot_temp = pot_temp + delta
342         endif
343        
# Line 394 | Line 354 | contains
354      fy = dudr * drdy
355      fz = dudr * drdz
356  
397 –
357   #ifdef IS_MPI
358      if (do_pot) then
359 <       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
360 <       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
359 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
360 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
361      endif
362  
363      f_Row(1,atom1) = f_Row(1,atom1) + fx

Diff Legend

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