ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/suttonchen.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/suttonchen.F90 (file contents):
Revision 2434 by chuckv, Tue Nov 15 16:18:36 2005 UTC vs.
Revision 2680 by chuckv, Thu Mar 30 21:08:48 2006 UTC

# Line 1 | Line 1
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
# Line 49 | Line 49 | module suttonchen
49    use status
50    use atype_module
51    use vector_class
52 +  use fForceOptions
53   #ifdef IS_MPI
54    use mpiSimulation
55   #endif
# Line 65 | Line 66 | module suttonchen
66    logical, save :: haveRcut = .false.
67    logical, save :: haveMixingMap = .false.
68    logical, save :: useGeometricDistanceMixing = .false.
69 +  logical, save :: cleanArrays = .true.
70 +  logical, save :: arraysAllocated = .false.
71  
72  
70 –
71 –
73    character(len = statusMsgSize) :: errMesg
74    integer :: sc_err
75  
76    character(len = 200) :: errMsg
77    character(len=*), parameter :: RoutineName =  "Sutton-Chen MODULE"
78    !! Logical that determines if eam arrays should be zeroed
78 –  logical :: cleanme = .true.
79    logical :: nmflag  = .false.
80  
81  
# Line 146 | Line 146 | module suttonchen
146    public :: getSCCut
147   ! public :: setSCDefaultCutoff
148   ! public :: setSCUniformCutoff
149 <  public :: useGeometricMixing
149 >
150  
151   contains
152  
# Line 176 | Line 176 | contains
176  
177      ! check to see if this is the first time into
178      if (.not.associated(SCList%SCTypes)) then
179 <       call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList)
179 >       call getMatchingElementList(atypes, "is_SC", .true., nSCtypes, MatchList)
180         SCList%nSCtypes = nSCtypes
181         allocate(SCList%SCTypes(nSCTypes))
182         nAtypes = getSize(atypes)
# Line 209 | Line 209 | contains
209         deallocate(SCList%atidToSCtype)
210         SCList%atidToSCtype=>null()
211      end if
212 + ! Reset Capacity
213 +    SCList% nSCTypes = 0
214 +    SCList%currentSCtype=0
215  
213 –
216    end subroutine destroySCTypes
217  
218  
# Line 221 | Line 223 | contains
223      real(kind=dp) :: cutValue
224      
225      scID = SCList%atidToSCType(atomID)
226 <    cutValue = SCList%SCTypes(scID)%sc_rcut
226 >    cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
227    end function getSCCut
228  
229  
# Line 243 | Line 245 | contains
245      if (.not. allocated(MixingMap)) then
246         allocate(MixingMap(nSCtypes, nSCtypes))
247      endif
248 <
248 >    useGeometricDistanceMixing = usesGeometricDistanceMixing()
249      do i = 1, nSCtypes
250  
251         e1 = SCList%SCtypes(i)%epsilon
# Line 272 | Line 274 | contains
274            MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
275            MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
276                 (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
277 +
278            if (i.ne.j) then
279               MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
280               MixingMap(j,i)%m          = MixingMap(i,j)%m
# Line 290 | Line 293 | contains
293  
294    !! routine checks to see if array is allocated, deallocates array if allocated
295    !! and then creates the array to the required size
296 <  subroutine allocateSC(status)
297 <    integer, intent(out) :: status
296 >  subroutine allocateSC()
297 >    integer :: status
298  
299   #ifdef IS_MPI
300      integer :: nAtomsInRow
# Line 299 | Line 302 | contains
302   #endif
303      integer :: alloc_stat
304  
305 <
305 >    
306      status = 0
307   #ifdef IS_MPI
308      nAtomsInRow = getNatomsInRow(plan_atom_row)
# Line 312 | Line 315 | contains
315      allocate(frho(nlocal),stat=alloc_stat)
316      if (alloc_stat /= 0) then
317         status = -1
315 –       return
318      end if
319  
320      if (allocated(rho)) deallocate(rho)
321      allocate(rho(nlocal),stat=alloc_stat)
322      if (alloc_stat /= 0) then
323         status = -1
322 –       return
324      end if
325  
326      if (allocated(dfrhodrho)) deallocate(dfrhodrho)
327      allocate(dfrhodrho(nlocal),stat=alloc_stat)
328      if (alloc_stat /= 0) then
329         status = -1
329 –       return
330      end if
331  
332   #ifdef IS_MPI
# Line 335 | Line 335 | contains
335      allocate(rho_tmp(nlocal),stat=alloc_stat)
336      if (alloc_stat /= 0) then
337         status = -1
338 –       return
338      end if
339  
340  
# Line 343 | Line 342 | contains
342      allocate(frho_row(nAtomsInRow),stat=alloc_stat)
343      if (alloc_stat /= 0) then
344         status = -1
346 –       return
345      end if
346      if (allocated(rho_row)) deallocate(rho_row)
347      allocate(rho_row(nAtomsInRow),stat=alloc_stat)
348      if (alloc_stat /= 0) then
349         status = -1
352 –       return
350      end if
351      if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
352      allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat)
353      if (alloc_stat /= 0) then
354         status = -1
358 –       return
355      end if
356  
357  
# Line 365 | Line 361 | contains
361      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
362      if (alloc_stat /= 0) then
363         status = -1
368 –       return
364      end if
365      if (allocated(rho_col)) deallocate(rho_col)
366      allocate(rho_col(nAtomsInCol),stat=alloc_stat)
367      if (alloc_stat /= 0) then
368         status = -1
374 –       return
369      end if
370      if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
371      allocate(dfrhodrho_col(nAtomsInCol),stat=alloc_stat)
372      if (alloc_stat /= 0) then
373         status = -1
380 –       return
374      end if
375  
376   #endif
377 <
377 >    if (status == -1) then
378 >       call handleError("SuttonChen:allocateSC","Error in allocating SC arrays")
379 >    end if
380 >    arraysAllocated = .true.
381    end subroutine allocateSC
382  
383    !! C sets rcut to be the largest cutoff of any atype
# Line 396 | Line 392 | contains
392  
393    end subroutine setCutoffSC
394  
395 <  subroutine useGeometricMixing()
400 <    useGeometricDistanceMixing = .true.
401 <    haveMixingMap = .false.
402 <    return
403 <  end subroutine useGeometricMixing
404 <  
405 <
406 <
407 <
408 <
409 <
410 <
411 <
412 <
395 > !! This array allocates module arrays if needed and builds mixing map.
396    subroutine clean_SC()
397 <
397 >    if (.not.arraysAllocated) call allocateSC()
398 >    if (.not.haveMixingMap) call createMixingMap()
399      ! clean non-IS_MPI first
400      frho = 0.0_dp
401      rho  = 0.0_dp
# Line 431 | Line 415 | contains
415  
416  
417    !! Calculates rho_r
418 <  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq)
418 >  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq, rcut)
419      integer :: atom1,atom2
420      real(kind = dp), dimension(3) :: d
421      real(kind = dp), intent(inout)               :: r
422 <    real(kind = dp), intent(inout)               :: rijsq
422 >    real(kind = dp), intent(inout)               :: rijsq, rcut
423      ! value of electron density rho do to atom i at atom j
424      real(kind = dp) :: rho_i_at_j
425      ! value of electron density rho do to atom j at atom i
# Line 452 | Line 436 | contains
436  
437      ! check to see if we need to be cleaned at the start of a force loop
438  
439 +    if (cleanArrays) call clean_SC()
440 +    cleanArrays = .false.
441  
456 –
457 –
442   #ifdef IS_MPI
443      Atid1 = Atid_row(Atom1)
444      Atid2 = Atid_col(Atom2)
# Line 489 | Line 473 | contains
473      real(kind=dp) :: pot
474      integer :: i,j
475      integer :: atom
492 –    real(kind=dp) :: U,U1,U2
476      integer :: atype1
477      integer :: atid1
478      integer :: myid
479  
480  
498 –    cleanme = .true.
481      !! Scatter the electron density from  pre-pair calculation back to local atoms
482   #ifdef IS_MPI
483      call scatter(rho_row,rho,plan_atom_row,sc_err)
# Line 519 | Line 501 | contains
501      do atom = 1, nlocal
502         Myid = SCList%atidtoSctype(Atid(atom))
503         frho(atom) = -SCList%SCTypes(Myid)%c * &
504 <            SCList%SCTypes(Myid)%epsilon * sqrt(rho(i))
504 >            SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom))
505  
506         dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
507 <       pot = pot + u
507 >
508 >       pot = pot + frho(atom)
509      enddo
510  
511   #ifdef IS_MPI
# Line 550 | Line 533 | contains
533  
534  
535    !! Does Sutton-Chen  pairwise Force calculation.  
536 <  subroutine do_sc_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
536 >  subroutine do_sc_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
537         pot, f, do_pot)
538      !Arguments    
539      integer, intent(in) ::  atom1, atom2
540 <    real( kind = dp ), intent(in) :: rij, r2
540 >    real( kind = dp ), intent(in) :: rij, r2, rcut
541      real( kind = dp ) :: pot, sw, vpair
542      real( kind = dp ), dimension(3,nLocal) :: f
543      real( kind = dp ), intent(in), dimension(3) :: d
# Line 566 | Line 549 | contains
549      real( kind = dp ) :: dvpdr
550      real( kind = dp ) :: drhodr
551      real( kind = dp ) :: dudr
569 –    real( kind = dp ) :: rcij
552      real( kind = dp ) :: drhoidr,drhojdr
553      real( kind = dp ) :: Fx,Fy,Fz
554 <    real( kind = dp ) :: r,d2pha,phb,d2phb
554 >    real( kind = dp ) :: d2pha,phb,d2phb
555      real( kind = dp ) :: pot_temp,vptmp
556      real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
557      integer :: id1,id2
# Line 579 | Line 561 | contains
561      !Local Variables
562  
563      ! write(*,*) "Frho: ", Frho(atom1)
564 +    
565 +    cleanArrays = .true.
566  
583 –
567      dvpdr = 0.0E0_DP
568  
569  
# Line 606 | Line 589 | contains
589         mij       = MixingMap(mytype_atom1,mytype_atom2)%m
590         vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
591  
592 <       vptmp = epsilonij*((aij/r)**nij)
592 >       vptmp = epsilonij*((aij/rij)**nij)
593  
594  
595 <       dvpdr = -nij*vptmp/r
596 <       drhodr = -mij*((aij/r)**mij)/r
595 >       dvpdr = -nij*vptmp/rij
596 >       drhodr = -mij*((aij/rij)**mij)/rij
597  
598        
599         dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &

Diff Legend

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