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 2534 by chuckv, Sat Dec 31 22:42:01 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

# Line 44 | Line 44 | module suttonchen
44  
45  
46   module suttonchen
47 +  use definitions
48    use simulation
49    use force_globals
50    use status
51    use atype_module
52    use vector_class
53    use fForceOptions
54 +  use interpolation
55   #ifdef IS_MPI
56    use mpiSimulation
57   #endif
# Line 58 | Line 60 | module suttonchen
60   #define __FORTRAN90
61   #include "UseTheForce/DarkSide/fInteractionMap.h"
62  
63 <  INTEGER, PARAMETER :: DP = selected_real_kind(15)
63 >  !! number of points for the spline approximations
64 >  INTEGER, PARAMETER :: np = 3000
65  
66    logical, save :: SC_FF_initialized = .false.
67    integer, save :: SC_Mixing_Policy
# Line 75 | Line 78 | module suttonchen
78  
79    character(len = 200) :: errMsg
80    character(len=*), parameter :: RoutineName =  "Sutton-Chen MODULE"
81 <  !! Logical that determines if eam arrays should be zeroed
79 <  logical :: nmflag  = .false.
80 <
81 <
81 >  
82    type, private :: SCtype
83 <     integer           :: atid      
84 <     real(kind=dp)     :: c
85 <     real(kind=dp)     :: m
86 <     real(kind=dp)     :: n
87 <     real(kind=dp)     :: alpha
88 <     real(kind=dp)     :: epsilon
89 <     real(kind=dp)     :: sc_rcut
83 >     integer       :: atid      
84 >     real(kind=dp) :: c
85 >     real(kind=dp) :: m
86 >     real(kind=dp) :: n
87 >     real(kind=dp) :: alpha
88 >     real(kind=dp) :: epsilon
89 >     real(kind=dp) :: sc_rcut
90    end type SCtype
91 +  
92  
92 –
93    !! Arrays for derivatives used in force calculation
94    real( kind = dp), dimension(:), allocatable :: rho
95    real( kind = dp), dimension(:), allocatable :: frho
96    real( kind = dp), dimension(:), allocatable :: dfrhodrho
97 <
98 <
99 <
97 >  
98    !! Arrays for MPI storage
99   #ifdef IS_MPI
100    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
# Line 110 | Line 108 | module suttonchen
108  
109    type, private :: SCTypeList
110       integer           :: nSCTypes = 0
111 <     integer           :: currentSCtype = 0
112 <
113 <     type (SCtype), pointer  :: SCtypes(:) => null()
116 <     integer, pointer         :: atidToSCtype(:) => null()
111 >     integer           :: currentSCtype = 0    
112 >     type (SCtype), pointer :: SCtypes(:) => null()
113 >     integer, pointer       :: atidToSCtype(:) => null()
114    end type SCTypeList
115  
116    type (SCTypeList), save :: SCList
117  
118 <
122 <
123 <
124 < type :: MixParameters
118 >  type:: MixParameters
119       real(kind=DP) :: alpha
120       real(kind=DP) :: epsilon
121 <     real(kind=DP) :: m
121 >     real(kind=DP) :: m
122       real(Kind=DP) :: n
129 –     real(kind=DP) :: vpair_pot
123       real(kind=dp) :: rCut
124 +     real(kind=dp) :: vCut
125       logical       :: rCutWasSet = .false.
126 <
126 >     type(cubicSpline) :: V
127 >     type(cubicSpline) :: phi
128    end type MixParameters
129  
130    type(MixParameters), dimension(:,:), allocatable :: MixingMap
131  
137 –
138 –
132    public :: setCutoffSC
133    public :: do_SC_pair
134    public :: newSCtype
# Line 157 | Line 150 | contains
150      real (kind = dp )                      :: n ! Pair Potential Exponent
151      real (kind = dp )                      :: alpha ! Length Scaling
152      real (kind = dp )                      :: epsilon ! Energy Scaling
160 –
161 –
153      integer                                :: c_ident
154      integer                                :: status
164 –
155      integer                                :: nAtypes,nSCTypes,myATID
156      integer                                :: maxVals
157      integer                                :: alloc_stat
# Line 188 | Line 178 | contains
178  
179      myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
180      SCList%atidToSCType(myATID) = current
191 –
192 –
181    
182      SCList%SCTypes(current)%atid         = c_ident
183      SCList%SCTypes(current)%alpha        = alpha
# Line 209 | Line 197 | contains
197         deallocate(SCList%atidToSCtype)
198         SCList%atidToSCtype=>null()
199      end if
200 + ! Reset Capacity
201 +    SCList%nSCTypes = 0
202 +    SCList%currentSCtype=0
203  
213 –
204    end subroutine destroySCTypes
205  
216 –
217 –
206    function getSCCut(atomID) result(cutValue)
207      integer, intent(in) :: atomID
208      integer :: scID
# Line 224 | Line 212 | contains
212      cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
213    end function getSCCut
214  
227 –
228 –
229 –
215    subroutine createMixingMap()
216 <    integer :: nSCtypes, i, j
217 <    real ( kind = dp ) :: e1, e2,m1,m2,alpha1,alpha2,n1,n2
218 <    real ( kind = dp ) :: rcut6, tp6, tp12
219 <    logical :: isSoftCore1, isSoftCore2, doShift
216 >    integer :: nSCtypes, i, j, k
217 >    real ( kind = dp ) :: e1, e2, m1, m2, alpha1, alpha2, n1, n2
218 >    real ( kind = dp ) :: epsilon, m, n, alpha, rCut, vCut, dr, r
219 >    real ( kind = dp ), dimension(np) :: rvals, vvals, phivals
220  
221      if (SCList%currentSCtype == 0) then
222         call handleError("SuttonChen", "No members in SCMap")
# Line 251 | Line 236 | contains
236         n1 = SCList%SCtypes(i)%n
237         alpha1 = SCList%SCtypes(i)%alpha
238  
239 <       do j = i, nSCtypes
239 >       do j = 1, nSCtypes
240            
256 –
241            e2 = SCList%SCtypes(j)%epsilon
242            m2 = SCList%SCtypes(j)%m
243            n2 = SCList%SCtypes(j)%n
244            alpha2 = SCList%SCtypes(j)%alpha
245  
246            if (useGeometricDistanceMixing) then
247 <             MixingMap(i,j)%alpha = sqrt(alpha1 * alpha2) !SC formulation
247 >             alpha = sqrt(alpha1 * alpha2) !SC formulation
248            else
249 <             MixingMap(i,j)%alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
249 >             alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
250            endif
251 +          rcut = 2.0_dp * alpha
252 +          epsilon = sqrt(e1 * e2)
253 +          m = 0.5_dp*(m1+m2)
254 +          n = 0.5_dp*(n1+n2)
255  
256 <          MixingMap(i,j)%epsilon = sqrt(e1 * e2)
257 <          MixingMap(i,j)%m = 0.5_dp*(m1+m2)
258 <          MixingMap(i,j)%n = 0.5_dp*(n1+n2)
259 <          MixingMap(i,j)%alpha = 0.5_dp*(alpha1+alpha2)
260 <          MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
261 <          MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
262 <               (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
263 <          if (i.ne.j) then
264 <             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
265 <             MixingMap(j,i)%m          = MixingMap(i,j)%m
266 <             MixingMap(j,i)%n          = MixingMap(i,j)%n
267 <             MixingMap(j,i)%alpha      = MixingMap(i,j)%alpha
268 <             MixingMap(j,i)%rcut = MixingMap(i,j)%rcut
269 <             MixingMap(j,i)%vpair_pot = MixingMap(i,j)%vpair_pot
270 <          endif
256 >          dr = (rCut) / dble(np-1)
257 >          rvals(1) = 0.0_dp
258 >          vvals(1) = 0.0_dp
259 >          phivals(1) = 0.0_dp
260 >
261 >          do k = 2, np
262 >             r = dble(k-1)*dr
263 >             rvals(k) = r
264 >             vvals(k) = epsilon*((alpha/r)**n)
265 >             phivals(k) = (alpha/r)**m
266 >          enddo
267 >
268 >          vCut = epsilon*((alpha/rCut)**n)
269 >
270 >          call newSpline(MixingMap(i,j)%V, rvals, vvals, .true.)
271 >          call newSpline(MixingMap(i,j)%phi, rvals, phivals, .true.)
272 >
273 >          MixingMap(i,j)%epsilon = epsilon
274 >          MixingMap(i,j)%m = m
275 >          MixingMap(i,j)%n = n
276 >          MixingMap(i,j)%alpha = alpha
277 >          MixingMap(i,j)%rCut = rcut
278 >          MixingMap(i,j)%vCut = vCut
279         enddo
280      enddo
281      
# Line 306 | Line 302 | contains
302      nAtomsInCol = getNatomsInCol(plan_atom_col)
303   #endif
304  
309 –
310 –
305      if (allocated(frho)) deallocate(frho)
306      allocate(frho(nlocal),stat=alloc_stat)
307      if (alloc_stat /= 0) then
# Line 377 | Line 371 | contains
371      arraysAllocated = .true.
372    end subroutine allocateSC
373  
374 <  !! C sets rcut to be the largest cutoff of any atype
381 <  !! present in this simulation. Doesn't include all atypes
382 <  !! sim knows about, just those in the simulation.
383 <  subroutine setCutoffSC(rcut, status)
374 >  subroutine setCutoffSC(rcut)
375      real(kind=dp) :: rcut
385 –    integer :: status
386 –    status = 0
387 –
376      SC_rcut = rcut
389 –
377    end subroutine setCutoffSC
378 <
379 < !! This array allocates module arrays if needed and builds mixing map.
378 >  
379 >  !! This array allocates module arrays if needed and builds mixing map.
380    subroutine clean_SC()
381      if (.not.arraysAllocated) call allocateSC()
382      if (.not.haveMixingMap) call createMixingMap()
# Line 409 | Line 396 | contains
396   #endif
397    end subroutine clean_SC
398  
412 –
413 –
399    !! Calculates rho_r
400 <  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq, rcut)
400 >  subroutine calc_sc_prepair_rho(atom1, atom2, d, r, rijsq, rcut)
401      integer :: atom1,atom2
402      real(kind = dp), dimension(3) :: d
403      real(kind = dp), intent(inout)               :: r
# Line 446 | Line 431 | contains
431  
432      Myid_atom1 = SCList%atidtoSCtype(Atid1)
433      Myid_atom2 = SCList%atidtoSCtype(Atid2)
434 +    
435 +    call lookupUniformSpline(MixingMap(myid_atom1,myid_atom2)%phi, r, &
436 +         rho_i_at_j)
437 +    rho_j_at_i = rho_i_at_j
438  
450 –
451 –
452 –       rho_i_at_j = (MixingMap(Myid_atom1,Myid_atom2)%alpha/r)&
453 –            **MixingMap(Myid_atom1,Myid_atom2)%m
454 –       rho_j_at_i = rho_i_at_j
455 –
439   #ifdef  IS_MPI
440 <       rho_col(atom2) = rho_col(atom2) + rho_i_at_j
441 <       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
440 >    rho_col(atom2) = rho_col(atom2) + rho_i_at_j
441 >    rho_row(atom1) = rho_row(atom1) + rho_j_at_i
442   #else
443 <       rho(atom2) = rho(atom2) + rho_i_at_j
444 <       rho(atom1) = rho(atom1) + rho_j_at_i
443 >    rho(atom2) = rho(atom2) + rho_i_at_j
444 >    rho(atom1) = rho(atom1) + rho_j_at_i
445   #endif
446 <
446 >    
447    end subroutine calc_sc_prepair_rho
448  
449  
450 < !! Calculate the rho_a for all local atoms
450 >  !! Calculate the rho_a for all local atoms
451    subroutine calc_sc_preforce_Frho(nlocal,pot)
452      integer :: nlocal
453      real(kind=dp) :: pot
454      integer :: i,j
455      integer :: atom
473 –    real(kind=dp) :: U,U1,U2
456      integer :: atype1
457      integer :: atid1
458      integer :: myid
459  
460 <
461 <    !! Scatter the electron density from  pre-pair calculation back to local atoms
460 >    !! Scatter the electron density from  pre-pair calculation back to
461 >    !! local atoms
462   #ifdef IS_MPI
463      call scatter(rho_row,rho,plan_atom_row,sc_err)
464      if (sc_err /= 0 ) then
# Line 492 | Line 474 | contains
474  
475      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
476   #endif
477 <
496 <
497 <
477 >    
478      !! Calculate F(rho) and derivative
479      do atom = 1, nlocal
480         Myid = SCList%atidtoSctype(Atid(atom))
481 <       frho(atom) = -SCList%SCTypes(Myid)%c * &
482 <            SCList%SCTypes(Myid)%epsilon * sqrt(rho(i))
481 >       frho(atom) = - SCList%SCTypes(Myid)%c * &
482 >            SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom))
483  
484         dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
485 <       pot = pot + u
485 >
486 >       pot = pot + frho(atom)
487      enddo
488  
489   #ifdef IS_MPI
# Line 526 | Line 507 | contains
507   #endif
508      
509      
510 <  end subroutine calc_sc_preforce_Frho
511 <
531 <
510 >  end subroutine calc_sc_preforce_Frho  
511 >  
512    !! Does Sutton-Chen  pairwise Force calculation.  
513    subroutine do_sc_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
514         pot, f, do_pot)
# Line 542 | Line 522 | contains
522  
523      logical, intent(in) :: do_pot
524  
525 <    real( kind = dp ) :: drdx,drdy,drdz
525 >    real( kind = dp ) :: drdx, drdy, drdz
526      real( kind = dp ) :: dvpdr
527 <    real( kind = dp ) :: drhodr
527 >    real( kind = dp ) :: rhtmp, drhodr
528      real( kind = dp ) :: dudr
549 –    real( kind = dp ) :: rcij
550 –    real( kind = dp ) :: drhoidr,drhojdr
529      real( kind = dp ) :: Fx,Fy,Fz
530 <    real( kind = dp ) :: r,d2pha,phb,d2phb
531 <    real( kind = dp ) :: pot_temp,vptmp
532 <    real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
555 <    integer :: id1,id2
530 >    real( kind = dp ) :: pot_temp, vptmp
531 >    real( kind = dp ) :: rcij, vcij
532 >    integer :: id1, id2
533      integer  :: mytype_atom1
534      integer  :: mytype_atom2
535 <    integer  :: atid1,atid2
535 >    integer  :: atid1, atid2
536      !Local Variables
560 –
561 –    ! write(*,*) "Frho: ", Frho(atom1)
537      
538      cleanArrays = .true.
539  
565 –    dvpdr = 0.0E0_DP
566 –
567 –
540   #ifdef IS_MPI
541         atid1 = atid_row(atom1)
542         atid2 = atid_col(atom2)
# Line 579 | Line 551 | contains
551         drdx = d(1)/rij
552         drdy = d(2)/rij
553         drdz = d(3)/rij
582 –
554                  
555 <       epsilonij = MixingMap(mytype_atom1,mytype_atom2)%epsilon
556 <       aij       = MixingMap(mytype_atom1,mytype_atom2)%alpha
586 <       nij       = MixingMap(mytype_atom1,mytype_atom2)%n
587 <       mij       = MixingMap(mytype_atom1,mytype_atom2)%m
588 <       vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
589 <
590 <       vptmp = epsilonij*((aij/r)**nij)
591 <
592 <
593 <       dvpdr = -nij*vptmp/r
594 <       drhodr = -mij*((aij/r)**mij)/r
595 <
555 >       rcij = MixingMap(mytype_atom1,mytype_atom2)%rCut
556 >       vcij = MixingMap(mytype_atom1,mytype_atom2)%vCut
557        
558 <       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
559 <            + dvpdr
599 <      
600 <       pot_temp = vptmp + vcij
601 <
558 >       call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%phi, &
559 >            rij, rhtmp, drhodr)
560  
561 +       call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%V, &
562 +            rij, vptmp, dvpdr)
563 +      
564   #ifdef IS_MPI
565 <       dudr = drhodr*(dfrhodrho_row(atom1)+dfrhodrho_col(atom2)) &
605 <            + dvpdr
606 <
565 >       dudr = drhodr*(dfrhodrho_row(atom1) + dfrhodrho_col(atom2)) + dvpdr
566   #else
567 <       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
609 <            + dvpdr
567 >       dudr = drhodr*(dfrhodrho(atom1)+ dfrhodrho(atom2)) + dvpdr
568   #endif
569 <
570 <
569 >              
570 >       pot_temp = vptmp - vcij
571 >
572         fx = dudr * drdx
573         fy = dudr * drdy
574         fz = dudr * drdz
575  
617 –
576   #ifdef IS_MPI
577         if (do_pot) then
578            pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + (pot_temp)*0.5
# Line 659 | Line 617 | contains
617            fpair(3) = fpair(3) + fz
618  
619         endif
662 –
663 –
620    end subroutine do_sc_pair
665 –
666 –
667 –
621   end module suttonchen

Diff Legend

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