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:	1388
Committed:	Fri Oct 30 16:38:48 2009 UTC (15 years, 8 months ago) by chuckv
Original Path:	*trunk/src/UseTheForce/DarkSide/suttonchen.F90*
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	939	!!
2	chuckv	1388	!! Copyright (c) 2005, 2009 The University of Notre Dame. All Rights Reserved.
3	chuckv	702	!!
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	960	use definitions
48	chuckv	702	use simulation
49			use force_globals
50			use status
51			use atype_module
52			use vector_class
53	chuckv	798	use fForceOptions
54	gezelter	938	use interpolation
55	chuckv	702	implicit none
56			PRIVATE
57			#define __FORTRAN90
58			#include "UseTheForce/DarkSide/fInteractionMap.h"
59
60	gezelter	938	!! number of points for the spline approximations
61			INTEGER, PARAMETER :: np = 3000
62	chuckv	702
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	728	logical, save :: haveMixingMap = .false.
68			logical, save :: useGeometricDistanceMixing = .false.
69	chuckv	835	logical, save :: cleanArrays = .true.
70			logical, save :: arraysAllocated = .false.
71	chuckv	1175
72	chuckv	702
73			character(len = statusMsgSize) :: errMesg
74	chuckv	714	integer :: sc_err
75	chuckv	702
76			character(len = 200) :: errMsg
77			character(len=*), parameter :: RoutineName = "Sutton-Chen MODULE"
78	gezelter	938
79	chuckv	702	type, private :: SCtype
80	gezelter	938	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	702	end type SCtype
88	gezelter	938
89	chuckv	702
90			type, private :: SCTypeList
91			integer :: nSCTypes = 0
92	gezelter	938	integer :: currentSCtype = 0
93			type (SCtype), pointer :: SCtypes(:) => null()
94			integer, pointer :: atidToSCtype(:) => null()
95	chuckv	702	end type SCTypeList
96
97			type (SCTypeList), save :: SCList
98
99	gezelter	938	type:: MixParameters
100	chuckv	707	real(kind=DP) :: alpha
101			real(kind=DP) :: epsilon
102	gezelter	938	real(kind=DP) :: m
103	chuckv	728	real(Kind=DP) :: n
104	chuckv	707	real(kind=dp) :: rCut
105	gezelter	938	real(kind=dp) :: vCut
106	chuckv	707	logical :: rCutWasSet = .false.
107	gezelter	938	type(cubicSpline) :: V
108			type(cubicSpline) :: phi
109	chuckv	707	end type MixParameters
110
111			type(MixParameters), dimension(:,:), allocatable :: MixingMap
112
113	chuckv	702	public :: setCutoffSC
114			public :: do_SC_pair
115			public :: newSCtype
116			public :: calc_SC_prepair_rho
117	chuckv	735	public :: calc_SC_preforce_Frho
118	chuckv	702	public :: destroySCtypes
119			public :: getSCCut
120	chuckv	728	! public :: setSCDefaultCutoff
121			! public :: setSCUniformCutoff
122	chuckv	798
123	chuckv	702
124			contains
125
126
127	chuckv	728	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	702	integer :: c_ident
134			integer :: status
135	chuckv	713	integer :: nAtypes,nSCTypes,myATID
136	chuckv	702	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	728	if (.not.associated(SCList%SCTypes)) then
149	chuckv	831	call getMatchingElementList(atypes, "is_SC", .true., nSCtypes, MatchList)
150	chuckv	728	SCList%nSCtypes = nSCtypes
151			allocate(SCList%SCTypes(nSCTypes))
152	chuckv	702	nAtypes = getSize(atypes)
153	chuckv	728	allocate(SCList%atidToSCType(nAtypes))
154	chuckv	1175	SCList%atidToSCType = -1
155	chuckv	702	end if
156
157	chuckv	728	SCList%currentSCType = SCList%currentSCType + 1
158			current = SCList%currentSCType
159	chuckv	702
160			myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
161	chuckv	728	SCList%atidToSCType(myATID) = current
162	chuckv	1175
163	chuckv	702
164	chuckv	728	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	713	end subroutine newSCtype
171	chuckv	702
172	chuckv	713
173			subroutine destroySCTypes()
174	chuckv	728	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	926	! Reset Capacity
183	gezelter	938	SCList%nSCTypes = 0
184	chuckv	926	SCList%currentSCtype=0
185	chuckv	702
186	chuckv	713	end subroutine destroySCTypes
187	chuckv	702
188	chuckv	713	function getSCCut(atomID) result(cutValue)
189	chuckv	702	integer, intent(in) :: atomID
190	chuckv	728	integer :: scID
191	chuckv	702	real(kind=dp) :: cutValue
192
193	chuckv	728	scID = SCList%atidToSCType(atomID)
194	chuckv	831	cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
195	chuckv	713	end function getSCCut
196	chuckv	702
197	chuckv	713	subroutine createMixingMap()
198	gezelter	938	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	713
203			if (SCList%currentSCtype == 0) then
204			call handleError("SuttonChen", "No members in SCMap")
205	chuckv	702	return
206			end if
207
208	chuckv	713	nSCtypes = SCList%nSCtypes
209	chuckv	702
210	chuckv	713	if (.not. allocated(MixingMap)) then
211			allocate(MixingMap(nSCtypes, nSCtypes))
212			endif
213	chuckv	798	useGeometricDistanceMixing = usesGeometricDistanceMixing()
214	chuckv	713	do i = 1, nSCtypes
215	chuckv	702
216	chuckv	713	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	702
221	gezelter	938	do j = 1, nSCtypes
222	chuckv	713
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	702
228	chuckv	728	if (useGeometricDistanceMixing) then
229	gezelter	938	alpha = sqrt(alpha1 * alpha2) !SC formulation
230	chuckv	728	else
231	gezelter	938	alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
232	chuckv	728	endif
233	gezelter	938	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	728
238	gezelter	938	dr = (rCut) / dble(np-1)
239	gezelter	960	rvals(1) = 0.0_dp
240			vvals(1) = 0.0_dp
241			phivals(1) = 0.0_dp
242	chuckv	842
243	gezelter	938	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	939	call newSpline(MixingMap(i,j)%V, rvals, vvals, .true.)
253			call newSpline(MixingMap(i,j)%phi, rvals, phivals, .true.)
254	gezelter	938
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	713	enddo
262	chuckv	702	enddo
263	chuckv	713
264			haveMixingMap = .true.
265
266			end subroutine createMixingMap
267
268	chuckv	702
269
270	gezelter	938	subroutine setCutoffSC(rcut)
271	chuckv	702	real(kind=dp) :: rcut
272	chuckv	713	SC_rcut = rcut
273			end subroutine setCutoffSC
274	gezelter	938
275	chuckv	702
276			!! Calculates rho_r
277	chuckv	1388	subroutine calc_sc_prepair_rho(atom1, atom2, atid1, atid2, d, r, rijsq, rho_i_at_j, rho_j_at_i, rcut)
278	gezelter	1386	integer :: atom1, atom2, atid1, atid2
279	chuckv	702	real(kind = dp), dimension(3) :: d
280			real(kind = dp), intent(inout) :: r
281	gezelter	762	real(kind = dp), intent(inout) :: rijsq, rcut
282	chuckv	702	! value of electron density rho do to atom i at atom j
283	chuckv	1388	real(kind = dp), intent(inout) :: rho_i_at_j
284	chuckv	702	! value of electron density rho do to atom j at atom i
285	chuckv	1388	real(kind = dp), intent(inout) :: rho_j_at_i
286	chuckv	702	! we don't use the derivatives, dummy variables
287	chuckv	728	real( kind = dp) :: drho
288	chuckv	714	integer :: sc_err
289	chuckv	702
290	gezelter	1285	integer :: myid_atom1 ! SC atid
291	chuckv	702	integer :: myid_atom2
292
293			! check to see if we need to be cleaned at the start of a force loop
294
295	chuckv	1388	if (.not.haveMixingMap) call createMixingMap()
296			haveMixingMap = .true.
297	chuckv	1175
298	chuckv	728	Myid_atom1 = SCList%atidtoSCtype(Atid1)
299			Myid_atom2 = SCList%atidtoSCtype(Atid2)
300	gezelter	938
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	1388
305	chuckv	713	end subroutine calc_sc_prepair_rho
306	chuckv	702
307
308	gezelter	938	!! Calculate the rho_a for all local atoms
309	gezelter	1285	subroutine calc_sc_preforce_Frho(nlocal, pot, particle_pot)
310	chuckv	702	integer :: nlocal
311			real(kind=dp) :: pot
312	gezelter	1285	real(kind=dp), dimension(nlocal) :: particle_pot
313	chuckv	702	integer :: i,j
314			integer :: atom
315			integer :: atype1
316	chuckv	728	integer :: atid1
317			integer :: myid
318	chuckv	1175
319	chuckv	713	!! Calculate F(rho) and derivative
320	chuckv	702	do atom = 1, nlocal
321	chuckv	728	Myid = SCList%atidtoSctype(Atid(atom))
322	chuckv	1175	! 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	728
331	chuckv	1175	dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
332			end if
333	chuckv	842	pot = pot + frho(atom)
334	gezelter	1285	particle_pot(atom) = particle_pot(atom) + frho(atom)
335	chuckv	702	enddo
336
337	gezelter	938	end subroutine calc_sc_preforce_Frho
338
339	chuckv	713	!! Does Sutton-Chen pairwise Force calculation.
340	gezelter	1386	subroutine do_sc_pair(atom1, atom2, atid1, atid2, d, rij, r2, rcut, sw, vpair, &
341	chuckv	1388	fpair, pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, fshift_i, fshift_j, do_pot)
342	chuckv	702	!Arguments
343			integer, intent(in) :: atom1, atom2
344	gezelter	762	real( kind = dp ), intent(in) :: rij, r2, rcut
345	chuckv	702	real( kind = dp ) :: pot, sw, vpair
346	gezelter	1386	real( kind = dp ), dimension(3) :: f1
347	chuckv	702	real( kind = dp ), intent(in), dimension(3) :: d
348			real( kind = dp ), intent(inout), dimension(3) :: fpair
349	chuckv	1388	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	702	logical, intent(in) :: do_pot
353
354	gezelter	938	real( kind = dp ) :: drdx, drdy, drdz
355	chuckv	728	real( kind = dp ) :: dvpdr
356	gezelter	938	real( kind = dp ) :: rhtmp, drhodr
357	chuckv	702	real( kind = dp ) :: dudr
358			real( kind = dp ) :: Fx,Fy,Fz
359	gezelter	938	real( kind = dp ) :: pot_temp, vptmp
360			real( kind = dp ) :: rcij, vcij
361			integer :: id1, id2
362	chuckv	702	integer :: mytype_atom1
363			integer :: mytype_atom2
364	gezelter	938	integer :: atid1, atid2
365	chuckv	702	!Local Variables
366	chuckv	835
367	gezelter	1386
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	1388
384			dudr = drhodr*(dfrhodrho_i + dfrhodrho_j) + dvpdr
385
386			pot_temp = vptmp - vcij
387	gezelter	1386
388	chuckv	1388	vpair = vpair + pot_temp
389	gezelter	938
390	chuckv	1388	fx = dudr * drdx
391			fy = dudr * drdy
392			fz = dudr * drdz
393	chuckv	702
394	chuckv	1388	if (do_pot) then
395	gezelter	1309
396	chuckv	1388	! 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	1309
416
417	chuckv	1388	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	702
423	chuckv	728	end subroutine do_sc_pair
424	chuckv	713	end module suttonchen