UseTheForce/DarkSide/suttonchen.F90

!!
!! Copyright (c) 2005, 2009 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Acknowledgement of the program authors must be made in any
!!    publication of scientific results based in part on use of the
!!    program.  An acceptable form of acknowledgement is citation of
!!    the article in which the program was described (Matthew
!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
!!    Parallel Simulation Engine for Molecular Dynamics,"
!!    J. Comput. Chem. 26, pp. 252-271 (2005))
!!
!! 2. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 3. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!

!! Impliments Sutton-Chen Metallic Potential
!! See A.P.SUTTON and J.CHEN,PHIL MAG LETT 61,139-146,1990


module suttonchen
  use definitions
  use simulation
  use force_globals
  use status
  use atype_module
  use vector_class
  use fForceOptions
  use interpolation
  implicit none
  PRIVATE
#define __FORTRAN90
#include "UseTheForce/DarkSide/fInteractionMap.h"

  !! number of points for the spline approximations
  INTEGER, PARAMETER :: np = 3000

  logical, save :: SC_FF_initialized = .false.
  integer, save :: SC_Mixing_Policy
  real(kind = dp), save :: SC_rcut
  logical, save :: haveRcut = .false.
  logical, save :: haveMixingMap = .false.
  logical, save :: useGeometricDistanceMixing = .false.
  logical, save :: cleanArrays = .true.
  logical, save :: arraysAllocated = .false.
  

  character(len = statusMsgSize) :: errMesg
  integer :: sc_err

  character(len = 200) :: errMsg
  character(len=*), parameter :: RoutineName =  "Sutton-Chen MODULE"
  
  type, private :: SCtype
     integer       :: atid       
     real(kind=dp) :: c 
     real(kind=dp) :: m
     real(kind=dp) :: n
     real(kind=dp) :: alpha
     real(kind=dp) :: epsilon
     real(kind=dp) :: sc_rcut
  end type SCtype
  

  type, private :: SCTypeList
     integer           :: nSCTypes = 0
     integer           :: currentSCtype = 0     
     type (SCtype), pointer :: SCtypes(:) => null()
     integer, pointer       :: atidToSCtype(:) => null()
  end type SCTypeList

  type (SCTypeList), save :: SCList

  type:: MixParameters
     real(kind=DP) :: alpha
     real(kind=DP) :: epsilon
     real(kind=DP) :: m 
     real(Kind=DP) :: n
     real(kind=dp) :: rCut
     real(kind=dp) :: vCut
     logical       :: rCutWasSet = .false.
     type(cubicSpline) :: V
     type(cubicSpline) :: phi
  end type MixParameters

  type(MixParameters), dimension(:,:), allocatable :: MixingMap

  public :: setCutoffSC
  public :: do_SC_pair
  public :: newSCtype
  public :: calc_SC_prepair_rho
  public :: calc_SC_preforce_Frho
  public :: destroySCtypes
  public :: getSCCut
 ! public :: setSCDefaultCutoff
 ! public :: setSCUniformCutoff
 

contains


  subroutine newSCtype(c_ident,c,m,n,alpha,epsilon,status)
    real (kind = dp )                      :: c ! Density Scaling
    real (kind = dp )                      :: m ! Density Exponent
    real (kind = dp )                      :: n ! Pair Potential Exponent
    real (kind = dp )                      :: alpha ! Length Scaling
    real (kind = dp )                      :: epsilon ! Energy Scaling
    integer                                :: c_ident
    integer                                :: status
    integer                                :: nAtypes,nSCTypes,myATID
    integer                                :: maxVals
    integer                                :: alloc_stat
    integer                                :: current
    integer,pointer                        :: Matchlist(:) => null()

    status = 0


    !! Assume that atypes has already been set and get the total number of types in atypes


    ! check to see if this is the first time into 
    if (.not.associated(SCList%SCTypes)) then
       call getMatchingElementList(atypes, "is_SC", .true., nSCtypes, MatchList)
       SCList%nSCtypes = nSCtypes
       allocate(SCList%SCTypes(nSCTypes))
       nAtypes = getSize(atypes)
       allocate(SCList%atidToSCType(nAtypes))
       SCList%atidToSCType = -1
    end if

    SCList%currentSCType = SCList%currentSCType + 1
    current = SCList%currentSCType

    myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
    SCList%atidToSCType(myATID) = current
    
  
    SCList%SCTypes(current)%atid         = c_ident
    SCList%SCTypes(current)%alpha        = alpha
    SCList%SCTypes(current)%c            = c
    SCList%SCTypes(current)%m            = m
    SCList%SCTypes(current)%n            = n
    SCList%SCTypes(current)%epsilon      = epsilon
  end subroutine newSCtype

  
  subroutine destroySCTypes()
    if (associated(SCList%SCtypes)) then
       deallocate(SCList%SCTypes)
       SCList%SCTypes=>null()
    end if
    if (associated(SCList%atidToSCtype)) then
       deallocate(SCList%atidToSCtype)
       SCList%atidToSCtype=>null()
    end if
! Reset Capacity
    SCList%nSCTypes = 0
    SCList%currentSCtype=0

  end subroutine destroySCTypes

  function getSCCut(atomID) result(cutValue)
    integer, intent(in) :: atomID
    integer :: scID
    real(kind=dp) :: cutValue
    
    scID = SCList%atidToSCType(atomID)
    cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
  end function getSCCut

  subroutine createMixingMap()
    integer :: nSCtypes, i, j, k
    real ( kind = dp ) :: e1, e2, m1, m2, alpha1, alpha2, n1, n2
    real ( kind = dp ) :: epsilon, m, n, alpha, rCut, vCut, dr, r
    real ( kind = dp ), dimension(np) :: rvals, vvals, phivals

    if (SCList%currentSCtype == 0) then
       call handleError("SuttonChen", "No members in SCMap")
       return
    end if

    nSCtypes = SCList%nSCtypes

    if (.not. allocated(MixingMap)) then
       allocate(MixingMap(nSCtypes, nSCtypes))
    endif
    useGeometricDistanceMixing = usesGeometricDistanceMixing()
    do i = 1, nSCtypes

       e1 = SCList%SCtypes(i)%epsilon
       m1 = SCList%SCtypes(i)%m
       n1 = SCList%SCtypes(i)%n
       alpha1 = SCList%SCtypes(i)%alpha

       do j = 1, nSCtypes
          
          e2 = SCList%SCtypes(j)%epsilon
          m2 = SCList%SCtypes(j)%m
          n2 = SCList%SCtypes(j)%n
          alpha2 = SCList%SCtypes(j)%alpha

          if (useGeometricDistanceMixing) then
             alpha = sqrt(alpha1 * alpha2) !SC formulation
          else
             alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
          endif
          rcut = 2.0_dp * alpha
          epsilon = sqrt(e1 * e2)
          m = 0.5_dp*(m1+m2)
          n = 0.5_dp*(n1+n2)

          dr = (rCut) / dble(np-1)
          rvals(1) = 0.0_dp
          vvals(1) = 0.0_dp
          phivals(1) = 0.0_dp

          do k = 2, np
             r = dble(k-1)*dr
             rvals(k) = r
             vvals(k) = epsilon*((alpha/r)**n)
             phivals(k) = (alpha/r)**m
          enddo

          vCut = epsilon*((alpha/rCut)**n)

          call newSpline(MixingMap(i,j)%V, rvals, vvals, .true.)
          call newSpline(MixingMap(i,j)%phi, rvals, phivals, .true.)

          MixingMap(i,j)%epsilon = epsilon
          MixingMap(i,j)%m = m
          MixingMap(i,j)%n = n
          MixingMap(i,j)%alpha = alpha
          MixingMap(i,j)%rCut = rcut
          MixingMap(i,j)%vCut = vCut
       enddo
    enddo
    
    haveMixingMap = .true.
    
  end subroutine createMixingMap
  

  subroutine setCutoffSC(rcut)
    real(kind=dp) :: rcut
    SC_rcut = rcut
  end subroutine setCutoffSC
  

  !! Calculates rho_r
  subroutine calc_sc_prepair_rho(atom1, atom2, atid1, atid2, d, r, rijsq, rho_i_at_j, rho_j_at_i, rcut)
    integer :: atom1, atom2, atid1, atid2
    real(kind = dp), dimension(3) :: d
    real(kind = dp), intent(inout)               :: r
    real(kind = dp), intent(inout)               :: rijsq, rcut
    ! value of electron density rho do to atom i at atom j
    real(kind = dp), intent(inout) :: rho_i_at_j
    ! value of electron density rho do to atom j at atom i
    real(kind = dp), intent(inout) :: rho_j_at_i
    ! we don't use the derivatives, dummy variables
    real( kind = dp) :: drho
    integer :: sc_err
    
    integer :: myid_atom1 ! SC atid
    integer :: myid_atom2 

    ! check to see if we need to be cleaned at the start of a force loop

    if (.not.haveMixingMap) call createMixingMap()
    haveMixingMap = .true.
    
    Myid_atom1 = SCList%atidtoSCtype(Atid1)
    Myid_atom2 = SCList%atidtoSCtype(Atid2)
    
    call lookupUniformSpline(MixingMap(myid_atom1,myid_atom2)%phi, r, &
         rho_i_at_j)
    rho_j_at_i = rho_i_at_j
       
  end subroutine calc_sc_prepair_rho


  !! Calculate the rho_a for all local atoms
  subroutine calc_sc_preforce_Frho(nlocal, pot, particle_pot)
    integer :: nlocal
    real(kind=dp) :: pot
    real(kind=dp), dimension(nlocal) :: particle_pot
    integer :: i,j
    integer :: atom
    integer :: atype1
    integer :: atid1
    integer :: myid
    
    !! Calculate F(rho) and derivative
    do atom = 1, nlocal
       Myid = SCList%atidtoSctype(Atid(atom))
       ! Myid is set to -1 for non SC atoms.
       ! Punt if we are a non-SC atom type.
       if (Myid == -1) then
          frho(atom) = 0.0_dp
          dfrhodrho(atom) = 0.0_dp
       else
          frho(atom) = - SCList%SCTypes(Myid)%c * &
              SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom))

          dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
       end if
       pot = pot + frho(atom)
       particle_pot(atom) = particle_pot(atom) + frho(atom)
    enddo

  end subroutine calc_sc_preforce_Frho  
  
  !! Does Sutton-Chen  pairwise Force calculation.  
  subroutine do_sc_pair(atom1, atom2, atid1, atid2, d, rij, r2, rcut, sw, vpair, &
        fpair, pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, fshift_i, fshift_j, do_pot)
    !Arguments    
    integer, intent(in) ::  atom1, atom2
    real( kind = dp ), intent(in) :: rij, r2, rcut
    real( kind = dp ) :: pot, sw, vpair
    real( kind = dp ), dimension(3) :: f1
    real( kind = dp ), intent(in), dimension(3) :: d
    real( kind = dp ), intent(inout), dimension(3) :: fpair
    real( kind = dp ), intent(inout) :: dfrhodrho_i, dfrhodrho_j 
    real( kind = dp ), intent(inout) :: rho_i, rho_j 
    real( kind = dp ), intent(inout):: fshift_i, fshift_j
    logical, intent(in) :: do_pot

    real( kind = dp ) :: drdx, drdy, drdz
    real( kind = dp ) :: dvpdr
    real( kind = dp ) :: rhtmp, drhodr
    real( kind = dp ) :: dudr
    real( kind = dp ) :: Fx,Fy,Fz
    real( kind = dp ) :: pot_temp, vptmp
    real( kind = dp ) :: rcij, vcij
    integer :: id1, id2
    integer  :: mytype_atom1
    integer  :: mytype_atom2
    integer  :: atid1, atid2
    !Local Variables
    
    
    mytype_atom1 = SCList%atidToSCType(atid1)
    mytype_atom2 = SCList%atidTOSCType(atid2)
    
    drdx = d(1)/rij
    drdy = d(2)/rij
    drdz = d(3)/rij
    
    rcij = MixingMap(mytype_atom1,mytype_atom2)%rCut
    vcij = MixingMap(mytype_atom1,mytype_atom2)%vCut
    
    call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%phi, &
         rij, rhtmp, drhodr)
    
    call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%V, &
         rij, vptmp, dvpdr)

    dudr = drhodr*(dfrhodrho_i + dfrhodrho_j) + dvpdr

    pot_temp = vptmp - vcij
    
    vpair = vpair + pot_temp
 
    fx = dudr * drdx
    fy = dudr * drdy
    fz = dudr * drdz

    if (do_pot) then

       ! particle_pot is the difference between the full potential 
       ! and the full potential without the presence of a particular
       ! particle (atom1).
       !
       ! This reduces the density at other particle locations, so
       ! we need to recompute the density at atom2 assuming atom1
       ! didn't contribute.  This then requires recomputing the
       ! density functional for atom2 as well.
       !
       ! Most of the particle_pot heavy lifting comes from the
       ! pair interaction, and will be handled by vpair.
       
       fshift_i = -SCList%SCTypes(mytype_atom1)%c * &
            SCList%SCTypes(mytype_atom1)%epsilon * &
            sqrt(rho_i-rhtmp)
       fshift_j = -SCList%SCTypes(mytype_atom2)%c * &
            SCList%SCTypes(mytype_atom2)%epsilon * &
            sqrt(rho_j-rhtmp)     
    end if


    pot = pot + pot_temp
    
    f1(1) = f1(1) + fx
    f1(2) = f1(2) + fy
    f1(3) = f1(3) + fz

  end subroutine do_sc_pair
end module suttonchen
Revision:	3561
Committed:	Fri Oct 30 16:38:48 2009 UTC (15 years, 10 months ago) by chuckv
File size:	13752 byte(s)
Log Message:	Removed remaining MPI from metallic potentials. Bug in MPI calculation of energies in Sutton-Chen.
#	User	Rev	Content
1	gezelter	2722	!!
2	chuckv	3561	!! Copyright (c) 2005, 2009 The University of Notre Dame. All Rights Reserved.
3	chuckv	2401	!!
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			!! Impliments Sutton-Chen Metallic Potential
43			!! See A.P.SUTTON and J.CHEN,PHIL MAG LETT 61,139-146,1990
44
45
46			module suttonchen
47	gezelter	2756	use definitions
48	chuckv	2401	use simulation
49			use force_globals
50			use status
51			use atype_module
52			use vector_class
53	chuckv	2497	use fForceOptions
54	gezelter	2717	use interpolation
55	chuckv	2401	implicit none
56			PRIVATE
57			#define __FORTRAN90
58			#include "UseTheForce/DarkSide/fInteractionMap.h"
59
60	gezelter	2717	!! number of points for the spline approximations
61			INTEGER, PARAMETER :: np = 3000
62	chuckv	2401
63			logical, save :: SC_FF_initialized = .false.
64			integer, save :: SC_Mixing_Policy
65			real(kind = dp), save :: SC_rcut
66			logical, save :: haveRcut = .false.
67	chuckv	2427	logical, save :: haveMixingMap = .false.
68			logical, save :: useGeometricDistanceMixing = .false.
69	chuckv	2534	logical, save :: cleanArrays = .true.
70			logical, save :: arraysAllocated = .false.
71	chuckv	3184
72	chuckv	2401
73			character(len = statusMsgSize) :: errMesg
74	chuckv	2413	integer :: sc_err
75	chuckv	2401
76			character(len = 200) :: errMsg
77			character(len=*), parameter :: RoutineName = "Sutton-Chen MODULE"
78	gezelter	2717
79	chuckv	2401	type, private :: SCtype
80	gezelter	2717	integer :: atid
81			real(kind=dp) :: c
82			real(kind=dp) :: m
83			real(kind=dp) :: n
84			real(kind=dp) :: alpha
85			real(kind=dp) :: epsilon
86			real(kind=dp) :: sc_rcut
87	chuckv	2401	end type SCtype
88	gezelter	2717
89	chuckv	2401
90			type, private :: SCTypeList
91			integer :: nSCTypes = 0
92	gezelter	2717	integer :: currentSCtype = 0
93			type (SCtype), pointer :: SCtypes(:) => null()
94			integer, pointer :: atidToSCtype(:) => null()
95	chuckv	2401	end type SCTypeList
96
97			type (SCTypeList), save :: SCList
98
99	gezelter	2717	type:: MixParameters
100	chuckv	2406	real(kind=DP) :: alpha
101			real(kind=DP) :: epsilon
102	gezelter	2717	real(kind=DP) :: m
103	chuckv	2427	real(Kind=DP) :: n
104	chuckv	2406	real(kind=dp) :: rCut
105	gezelter	2717	real(kind=dp) :: vCut
106	chuckv	2406	logical :: rCutWasSet = .false.
107	gezelter	2717	type(cubicSpline) :: V
108			type(cubicSpline) :: phi
109	chuckv	2406	end type MixParameters
110
111			type(MixParameters), dimension(:,:), allocatable :: MixingMap
112
113	chuckv	2401	public :: setCutoffSC
114			public :: do_SC_pair
115			public :: newSCtype
116			public :: calc_SC_prepair_rho
117	chuckv	2434	public :: calc_SC_preforce_Frho
118	chuckv	2401	public :: destroySCtypes
119			public :: getSCCut
120	chuckv	2427	! public :: setSCDefaultCutoff
121			! public :: setSCUniformCutoff
122	chuckv	2497
123	chuckv	2401
124			contains
125
126
127	chuckv	2427	subroutine newSCtype(c_ident,c,m,n,alpha,epsilon,status)
128			real (kind = dp ) :: c ! Density Scaling
129			real (kind = dp ) :: m ! Density Exponent
130			real (kind = dp ) :: n ! Pair Potential Exponent
131			real (kind = dp ) :: alpha ! Length Scaling
132			real (kind = dp ) :: epsilon ! Energy Scaling
133	chuckv	2401	integer :: c_ident
134			integer :: status
135	chuckv	2412	integer :: nAtypes,nSCTypes,myATID
136	chuckv	2401	integer :: maxVals
137			integer :: alloc_stat
138			integer :: current
139			integer,pointer :: Matchlist(:) => null()
140
141			status = 0
142
143
144			!! Assume that atypes has already been set and get the total number of types in atypes
145
146
147			! check to see if this is the first time into
148	chuckv	2427	if (.not.associated(SCList%SCTypes)) then
149	chuckv	2530	call getMatchingElementList(atypes, "is_SC", .true., nSCtypes, MatchList)
150	chuckv	2427	SCList%nSCtypes = nSCtypes
151			allocate(SCList%SCTypes(nSCTypes))
152	chuckv	2401	nAtypes = getSize(atypes)
153	chuckv	2427	allocate(SCList%atidToSCType(nAtypes))
154	chuckv	3184	SCList%atidToSCType = -1
155	chuckv	2401	end if
156
157	chuckv	2427	SCList%currentSCType = SCList%currentSCType + 1
158			current = SCList%currentSCType
159	chuckv	2401
160			myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
161	chuckv	2427	SCList%atidToSCType(myATID) = current
162	chuckv	3184
163	chuckv	2401
164	chuckv	2427	SCList%SCTypes(current)%atid = c_ident
165			SCList%SCTypes(current)%alpha = alpha
166			SCList%SCTypes(current)%c = c
167			SCList%SCTypes(current)%m = m
168			SCList%SCTypes(current)%n = n
169			SCList%SCTypes(current)%epsilon = epsilon
170	chuckv	2412	end subroutine newSCtype
171	chuckv	2401
172	chuckv	2412
173			subroutine destroySCTypes()
174	chuckv	2427	if (associated(SCList%SCtypes)) then
175			deallocate(SCList%SCTypes)
176			SCList%SCTypes=>null()
177			end if
178			if (associated(SCList%atidToSCtype)) then
179			deallocate(SCList%atidToSCtype)
180			SCList%atidToSCtype=>null()
181			end if
182	chuckv	2680	! Reset Capacity
183	gezelter	2717	SCList%nSCTypes = 0
184	chuckv	2680	SCList%currentSCtype=0
185	chuckv	2401
186	chuckv	2412	end subroutine destroySCTypes
187	chuckv	2401
188	chuckv	2412	function getSCCut(atomID) result(cutValue)
189	chuckv	2401	integer, intent(in) :: atomID
190	chuckv	2427	integer :: scID
191	chuckv	2401	real(kind=dp) :: cutValue
192
193	chuckv	2427	scID = SCList%atidToSCType(atomID)
194	chuckv	2530	cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
195	chuckv	2412	end function getSCCut
196	chuckv	2401
197	chuckv	2412	subroutine createMixingMap()
198	gezelter	2717	integer :: nSCtypes, i, j, k
199			real ( kind = dp ) :: e1, e2, m1, m2, alpha1, alpha2, n1, n2
200			real ( kind = dp ) :: epsilon, m, n, alpha, rCut, vCut, dr, r
201			real ( kind = dp ), dimension(np) :: rvals, vvals, phivals
202	chuckv	2412
203			if (SCList%currentSCtype == 0) then
204			call handleError("SuttonChen", "No members in SCMap")
205	chuckv	2401	return
206			end if
207
208	chuckv	2412	nSCtypes = SCList%nSCtypes
209	chuckv	2401
210	chuckv	2412	if (.not. allocated(MixingMap)) then
211			allocate(MixingMap(nSCtypes, nSCtypes))
212			endif
213	chuckv	2497	useGeometricDistanceMixing = usesGeometricDistanceMixing()
214	chuckv	2412	do i = 1, nSCtypes
215	chuckv	2401
216	chuckv	2412	e1 = SCList%SCtypes(i)%epsilon
217			m1 = SCList%SCtypes(i)%m
218			n1 = SCList%SCtypes(i)%n
219			alpha1 = SCList%SCtypes(i)%alpha
220	chuckv	2401
221	gezelter	2717	do j = 1, nSCtypes
222	chuckv	2412
223			e2 = SCList%SCtypes(j)%epsilon
224			m2 = SCList%SCtypes(j)%m
225			n2 = SCList%SCtypes(j)%n
226			alpha2 = SCList%SCtypes(j)%alpha
227	chuckv	2401
228	chuckv	2427	if (useGeometricDistanceMixing) then
229	gezelter	2717	alpha = sqrt(alpha1 * alpha2) !SC formulation
230	chuckv	2427	else
231	gezelter	2717	alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
232	chuckv	2427	endif
233	gezelter	2717	rcut = 2.0_dp * alpha
234			epsilon = sqrt(e1 * e2)
235			m = 0.5_dp*(m1+m2)
236			n = 0.5_dp*(n1+n2)
237	chuckv	2427
238	gezelter	2717	dr = (rCut) / dble(np-1)
239	gezelter	2756	rvals(1) = 0.0_dp
240			vvals(1) = 0.0_dp
241			phivals(1) = 0.0_dp
242	chuckv	2541
243	gezelter	2717	do k = 2, np
244			r = dble(k-1)*dr
245			rvals(k) = r
246			vvals(k) = epsilon((alpha/r)*n)
247			phivals(k) = (alpha/r)**m
248			enddo
249
250			vCut = epsilon((alpha/rCut)*n)
251
252	gezelter	2722	call newSpline(MixingMap(i,j)%V, rvals, vvals, .true.)
253			call newSpline(MixingMap(i,j)%phi, rvals, phivals, .true.)
254	gezelter	2717
255			MixingMap(i,j)%epsilon = epsilon
256			MixingMap(i,j)%m = m
257			MixingMap(i,j)%n = n
258			MixingMap(i,j)%alpha = alpha
259			MixingMap(i,j)%rCut = rcut
260			MixingMap(i,j)%vCut = vCut
261	chuckv	2412	enddo
262	chuckv	2401	enddo
263	chuckv	2412
264			haveMixingMap = .true.
265
266			end subroutine createMixingMap
267
268	chuckv	2401
269
270	gezelter	2717	subroutine setCutoffSC(rcut)
271	chuckv	2401	real(kind=dp) :: rcut
272	chuckv	2412	SC_rcut = rcut
273			end subroutine setCutoffSC
274	gezelter	2717
275	chuckv	2401
276			!! Calculates rho_r
277	chuckv	3561	subroutine calc_sc_prepair_rho(atom1, atom2, atid1, atid2, d, r, rijsq, rho_i_at_j, rho_j_at_i, rcut)
278	gezelter	3559	integer :: atom1, atom2, atid1, atid2
279	chuckv	2401	real(kind = dp), dimension(3) :: d
280			real(kind = dp), intent(inout) :: r
281	gezelter	2461	real(kind = dp), intent(inout) :: rijsq, rcut
282	chuckv	2401	! value of electron density rho do to atom i at atom j
283	chuckv	3561	real(kind = dp), intent(inout) :: rho_i_at_j
284	chuckv	2401	! value of electron density rho do to atom j at atom i
285	chuckv	3561	real(kind = dp), intent(inout) :: rho_j_at_i
286	chuckv	2401	! we don't use the derivatives, dummy variables
287	chuckv	2427	real( kind = dp) :: drho
288	chuckv	2413	integer :: sc_err
289	chuckv	2401
290	gezelter	3440	integer :: myid_atom1 ! SC atid
291	chuckv	2401	integer :: myid_atom2
292
293			! check to see if we need to be cleaned at the start of a force loop
294
295	chuckv	3561	if (.not.haveMixingMap) call createMixingMap()
296			haveMixingMap = .true.
297	chuckv	3184
298	chuckv	2427	Myid_atom1 = SCList%atidtoSCtype(Atid1)
299			Myid_atom2 = SCList%atidtoSCtype(Atid2)
300	gezelter	2717
301			call lookupUniformSpline(MixingMap(myid_atom1,myid_atom2)%phi, r, &
302			rho_i_at_j)
303			rho_j_at_i = rho_i_at_j
304	chuckv	3561
305	chuckv	2412	end subroutine calc_sc_prepair_rho
306	chuckv	2401
307
308	gezelter	2717	!! Calculate the rho_a for all local atoms
309	gezelter	3440	subroutine calc_sc_preforce_Frho(nlocal, pot, particle_pot)
310	chuckv	2401	integer :: nlocal
311			real(kind=dp) :: pot
312	gezelter	3440	real(kind=dp), dimension(nlocal) :: particle_pot
313	chuckv	2401	integer :: i,j
314			integer :: atom
315			integer :: atype1
316	chuckv	2427	integer :: atid1
317			integer :: myid
318	chuckv	3184
319	chuckv	2412	!! Calculate F(rho) and derivative
320	chuckv	2401	do atom = 1, nlocal
321	chuckv	2427	Myid = SCList%atidtoSctype(Atid(atom))
322	chuckv	3184	! Myid is set to -1 for non SC atoms.
323			! Punt if we are a non-SC atom type.
324			if (Myid == -1) then
325			frho(atom) = 0.0_dp
326			dfrhodrho(atom) = 0.0_dp
327			else
328			frho(atom) = - SCList%SCTypes(Myid)%c * &
329			SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom))
330	chuckv	2427
331	chuckv	3184	dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
332			end if
333	chuckv	2541	pot = pot + frho(atom)
334	gezelter	3440	particle_pot(atom) = particle_pot(atom) + frho(atom)
335	chuckv	2401	enddo
336
337	gezelter	2717	end subroutine calc_sc_preforce_Frho
338
339	chuckv	2412	!! Does Sutton-Chen pairwise Force calculation.
340	gezelter	3559	subroutine do_sc_pair(atom1, atom2, atid1, atid2, d, rij, r2, rcut, sw, vpair, &
341	chuckv	3561	fpair, pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, fshift_i, fshift_j, do_pot)
342	chuckv	2401	!Arguments
343			integer, intent(in) :: atom1, atom2
344	gezelter	2461	real( kind = dp ), intent(in) :: rij, r2, rcut
345	chuckv	2401	real( kind = dp ) :: pot, sw, vpair
346	gezelter	3559	real( kind = dp ), dimension(3) :: f1
347	chuckv	2401	real( kind = dp ), intent(in), dimension(3) :: d
348			real( kind = dp ), intent(inout), dimension(3) :: fpair
349	chuckv	3561	real( kind = dp ), intent(inout) :: dfrhodrho_i, dfrhodrho_j
350			real( kind = dp ), intent(inout) :: rho_i, rho_j
351			real( kind = dp ), intent(inout):: fshift_i, fshift_j
352	chuckv	2401	logical, intent(in) :: do_pot
353
354	gezelter	2717	real( kind = dp ) :: drdx, drdy, drdz
355	chuckv	2427	real( kind = dp ) :: dvpdr
356	gezelter	2717	real( kind = dp ) :: rhtmp, drhodr
357	chuckv	2401	real( kind = dp ) :: dudr
358			real( kind = dp ) :: Fx,Fy,Fz
359	gezelter	2717	real( kind = dp ) :: pot_temp, vptmp
360			real( kind = dp ) :: rcij, vcij
361			integer :: id1, id2
362	chuckv	2401	integer :: mytype_atom1
363			integer :: mytype_atom2
364	gezelter	2717	integer :: atid1, atid2
365	chuckv	2401	!Local Variables
366	chuckv	2534
367	gezelter	3559
368			mytype_atom1 = SCList%atidToSCType(atid1)
369			mytype_atom2 = SCList%atidTOSCType(atid2)
370
371			drdx = d(1)/rij
372			drdy = d(2)/rij
373			drdz = d(3)/rij
374
375			rcij = MixingMap(mytype_atom1,mytype_atom2)%rCut
376			vcij = MixingMap(mytype_atom1,mytype_atom2)%vCut
377
378			call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%phi, &
379			rij, rhtmp, drhodr)
380
381			call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%V, &
382			rij, vptmp, dvpdr)
383	chuckv	3561
384			dudr = drhodr*(dfrhodrho_i + dfrhodrho_j) + dvpdr
385
386			pot_temp = vptmp - vcij
387	gezelter	3559
388	chuckv	3561	vpair = vpair + pot_temp
389	gezelter	2717
390	chuckv	3561	fx = dudr * drdx
391			fy = dudr * drdy
392			fz = dudr * drdz
393	chuckv	2401
394	chuckv	3561	if (do_pot) then
395	gezelter	3466
396	chuckv	3561	! particle_pot is the difference between the full potential
397			! and the full potential without the presence of a particular
398			! particle (atom1).
399			!
400			! This reduces the density at other particle locations, so
401			! we need to recompute the density at atom2 assuming atom1
402			! didn't contribute. This then requires recomputing the
403			! density functional for atom2 as well.
404			!
405			! Most of the particle_pot heavy lifting comes from the
406			! pair interaction, and will be handled by vpair.
407
408			fshift_i = -SCList%SCTypes(mytype_atom1)%c * &
409			SCList%SCTypes(mytype_atom1)%epsilon * &
410			sqrt(rho_i-rhtmp)
411			fshift_j = -SCList%SCTypes(mytype_atom2)%c * &
412			SCList%SCTypes(mytype_atom2)%epsilon * &
413			sqrt(rho_j-rhtmp)
414			end if
415	gezelter	3466
416
417	chuckv	3561	pot = pot + pot_temp
418
419			f1(1) = f1(1) + fx
420			f1(2) = f1(2) + fy
421			f1(3) = f1(3) + fz
422	chuckv	2401
423	chuckv	2427	end subroutine do_sc_pair
424	chuckv	2412	end module suttonchen