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. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 2. 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.
!!
!! SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
!! research, please cite the appropriate papers when you publish your
!! work.  Good starting points are:
!!                                                                      
!! [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
!! [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
!! [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
!! [4]  Vardeman & Gezelter, in progress (2009).
!!

!! Implements 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(atid1, atid2, d, r, rijsq, rho_i_at_j, rho_j_at_i)
    integer :: atid1, atid2
    real(kind = dp), dimension(3) :: d
    real(kind = dp), intent(inout)               :: r
    real(kind = dp), intent(inout)               :: rijsq
    ! 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(atid, rho, frho, dfrhodrho)
    real(kind=dp) :: rho, frho, dfrhodrho
    real(kind=dp) :: u, u1
    integer :: atid, myid
    
    !! Calculate F(rho) and derivative
    
    Myid = SCList%atidtoSctype(atid)
    ! Myid is set to -1 for non SC atoms.
    ! Punt if we are a non-SC atom type.
    if (Myid == -1) then
       frho = 0.0_dp
       dfrhodrho = 0.0_dp
    else
       frho = - SCList%SCTypes(Myid)%c * &
            SCList%SCTypes(Myid)%epsilon * sqrt(rho)
          
       dfrhodrho = 0.5_dp*frho/rho
    end if
    
  end subroutine calc_sc_preforce_Frho
  
  !! Does Sutton-Chen  pairwise Force calculation.  
  subroutine do_sc_pair(atid1, atid2, d, rij, r2, sw, vpair, &
       pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, &
       fshift_i, fshift_j)
    !Arguments    
    integer, intent(in) ::  atid1, atid2
    real( kind = dp ), intent(in) :: rij, r2
    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) :: dfrhodrho_i, dfrhodrho_j 
    real( kind = dp ), intent(inout) :: rho_i, rho_j 
    real( kind = dp ), intent(inout):: fshift_i, fshift_j   
    
    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
    
    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
        
    ! 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)         
    
    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:	1479
Committed:	Mon Jul 26 19:00:48 2010 UTC (15 years ago) by gezelter
File size:	13097 byte(s)
Log Message:	Added EAM. Still segfaults but compiles.
#	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	gezelter	1390	!! 1. Redistributions of source code must retain the above copyright
10	chuckv	702	!! notice, this list of conditions and the following disclaimer.
11			!!
12	gezelter	1390	!! 2. Redistributions in binary form must reproduce the above copyright
13	chuckv	702	!! notice, this list of conditions and the following disclaimer in the
14			!! documentation and/or other materials provided with the
15			!! distribution.
16			!!
17			!! This software is provided "AS IS," without a warranty of any
18			!! kind. All express or implied conditions, representations and
19			!! warranties, including any implied warranty of merchantability,
20			!! fitness for a particular purpose or non-infringement, are hereby
21			!! excluded. The University of Notre Dame and its licensors shall not
22			!! be liable for any damages suffered by licensee as a result of
23			!! using, modifying or distributing the software or its
24			!! derivatives. In no event will the University of Notre Dame or its
25			!! licensors be liable for any lost revenue, profit or data, or for
26			!! direct, indirect, special, consequential, incidental or punitive
27			!! damages, however caused and regardless of the theory of liability,
28			!! arising out of the use of or inability to use software, even if the
29			!! University of Notre Dame has been advised of the possibility of
30			!! such damages.
31			!!
32	gezelter	1390	!! SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			!! research, please cite the appropriate papers when you publish your
34			!! work. Good starting points are:
35			!!
36			!! [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			!! [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			!! [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39			!! [4] Vardeman & Gezelter, in progress (2009).
40			!!
41	chuckv	702
42	gezelter	1448	!! Implements Sutton-Chen Metallic Potential
43	chuckv	702	!! 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	gezelter	1448
60	gezelter	938	!! number of points for the spline approximations
61			INTEGER, PARAMETER :: np = 3000
62	gezelter	1448
63	chuckv	702	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	gezelter	1448
76	chuckv	702	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	gezelter	1448
90	chuckv	702	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	gezelter	1448
97	chuckv	702	type (SCTypeList), save :: SCList
98	gezelter	1448
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	gezelter	1448
111	chuckv	707	type(MixParameters), dimension(:,:), allocatable :: MixingMap
112	gezelter	1448
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	gezelter	1448	! public :: setSCDefaultCutoff
121			! public :: setSCUniformCutoff
122
123
124	chuckv	702	contains
125	gezelter	1448
126
127	chuckv	728	subroutine newSCtype(c_ident,c,m,n,alpha,epsilon,status)
128	gezelter	1448	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			integer :: c_ident
134			integer :: status
135			integer :: nAtypes,nSCTypes,myATID
136			integer :: maxVals
137			integer :: alloc_stat
138			integer :: current
139			integer,pointer :: Matchlist(:) => null()
140
141	chuckv	702	status = 0
142	gezelter	1448
143
144	chuckv	702	!! Assume that atypes has already been set and get the total number of types in atypes
145	gezelter	1448
146	chuckv	702	! check to see if this is the first time into
147	chuckv	728	if (.not.associated(SCList%SCTypes)) then
148	chuckv	831	call getMatchingElementList(atypes, "is_SC", .true., nSCtypes, MatchList)
149	chuckv	728	SCList%nSCtypes = nSCtypes
150			allocate(SCList%SCTypes(nSCTypes))
151	chuckv	702	nAtypes = getSize(atypes)
152	chuckv	728	allocate(SCList%atidToSCType(nAtypes))
153	chuckv	1175	SCList%atidToSCType = -1
154	chuckv	702	end if
155	gezelter	1448
156	chuckv	728	SCList%currentSCType = SCList%currentSCType + 1
157			current = SCList%currentSCType
158	gezelter	1448
159	chuckv	702	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
160	chuckv	728	SCList%atidToSCType(myATID) = current
161	chuckv	1175
162	gezelter	1448
163	chuckv	728	SCList%SCTypes(current)%atid = c_ident
164			SCList%SCTypes(current)%alpha = alpha
165			SCList%SCTypes(current)%c = c
166			SCList%SCTypes(current)%m = m
167			SCList%SCTypes(current)%n = n
168			SCList%SCTypes(current)%epsilon = epsilon
169	chuckv	713	end subroutine newSCtype
170
171	gezelter	1448
172	chuckv	713	subroutine destroySCTypes()
173	chuckv	728	if (associated(SCList%SCtypes)) then
174			deallocate(SCList%SCTypes)
175			SCList%SCTypes=>null()
176			end if
177			if (associated(SCList%atidToSCtype)) then
178			deallocate(SCList%atidToSCtype)
179			SCList%atidToSCtype=>null()
180			end if
181	gezelter	1448	! Reset Capacity
182	gezelter	938	SCList%nSCTypes = 0
183	chuckv	926	SCList%currentSCtype=0
184	gezelter	1448
185	chuckv	713	end subroutine destroySCTypes
186	gezelter	1448
187	chuckv	713	function getSCCut(atomID) result(cutValue)
188	chuckv	702	integer, intent(in) :: atomID
189	chuckv	728	integer :: scID
190	chuckv	702	real(kind=dp) :: cutValue
191
192	chuckv	728	scID = SCList%atidToSCType(atomID)
193	chuckv	831	cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
194	chuckv	713	end function getSCCut
195	gezelter	1448
196	chuckv	713	subroutine createMixingMap()
197	gezelter	938	integer :: nSCtypes, i, j, k
198			real ( kind = dp ) :: e1, e2, m1, m2, alpha1, alpha2, n1, n2
199			real ( kind = dp ) :: epsilon, m, n, alpha, rCut, vCut, dr, r
200			real ( kind = dp ), dimension(np) :: rvals, vvals, phivals
201	gezelter	1448
202	chuckv	713	if (SCList%currentSCtype == 0) then
203			call handleError("SuttonChen", "No members in SCMap")
204	chuckv	702	return
205			end if
206	gezelter	1448
207	chuckv	713	nSCtypes = SCList%nSCtypes
208	gezelter	1448
209	chuckv	713	if (.not. allocated(MixingMap)) then
210			allocate(MixingMap(nSCtypes, nSCtypes))
211			endif
212	chuckv	798	useGeometricDistanceMixing = usesGeometricDistanceMixing()
213	chuckv	713	do i = 1, nSCtypes
214	gezelter	1448
215	chuckv	713	e1 = SCList%SCtypes(i)%epsilon
216			m1 = SCList%SCtypes(i)%m
217			n1 = SCList%SCtypes(i)%n
218			alpha1 = SCList%SCtypes(i)%alpha
219	gezelter	1448
220	gezelter	938	do j = 1, nSCtypes
221	chuckv	713
222			e2 = SCList%SCtypes(j)%epsilon
223			m2 = SCList%SCtypes(j)%m
224			n2 = SCList%SCtypes(j)%n
225			alpha2 = SCList%SCtypes(j)%alpha
226	gezelter	1448
227	chuckv	728	if (useGeometricDistanceMixing) then
228	gezelter	1448	alpha = sqrt(alpha1 * alpha2) ! SC formulation
229	chuckv	728	else
230	gezelter	938	alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
231	chuckv	728	endif
232	gezelter	938	rcut = 2.0_dp * alpha
233			epsilon = sqrt(e1 * e2)
234			m = 0.5_dp*(m1+m2)
235			n = 0.5_dp*(n1+n2)
236	gezelter	1448
237	gezelter	938	dr = (rCut) / dble(np-1)
238	gezelter	960	rvals(1) = 0.0_dp
239			vvals(1) = 0.0_dp
240			phivals(1) = 0.0_dp
241	gezelter	1448
242	gezelter	938	do k = 2, np
243			r = dble(k-1)*dr
244			rvals(k) = r
245			vvals(k) = epsilon((alpha/r)*n)
246			phivals(k) = (alpha/r)**m
247			enddo
248	gezelter	1448
249	gezelter	938	vCut = epsilon((alpha/rCut)*n)
250	gezelter	1448
251	gezelter	939	call newSpline(MixingMap(i,j)%V, rvals, vvals, .true.)
252			call newSpline(MixingMap(i,j)%phi, rvals, phivals, .true.)
253	gezelter	1448
254	gezelter	938	MixingMap(i,j)%epsilon = epsilon
255			MixingMap(i,j)%m = m
256			MixingMap(i,j)%n = n
257			MixingMap(i,j)%alpha = alpha
258			MixingMap(i,j)%rCut = rcut
259			MixingMap(i,j)%vCut = vCut
260	chuckv	713	enddo
261	chuckv	702	enddo
262	chuckv	713
263			haveMixingMap = .true.
264
265			end subroutine createMixingMap
266
267	chuckv	702
268	gezelter	1448
269	gezelter	938	subroutine setCutoffSC(rcut)
270	chuckv	702	real(kind=dp) :: rcut
271	chuckv	713	SC_rcut = rcut
272			end subroutine setCutoffSC
273	gezelter	938
274	gezelter	1448
275	chuckv	702	!! Calculates rho_r
276	gezelter	1464	subroutine calc_sc_prepair_rho(atid1, atid2, d, r, rijsq, rho_i_at_j, rho_j_at_i)
277			integer :: atid1, atid2
278	chuckv	702	real(kind = dp), dimension(3) :: d
279			real(kind = dp), intent(inout) :: r
280	gezelter	1419	real(kind = dp), intent(inout) :: rijsq
281	chuckv	702	! value of electron density rho do to atom i at atom j
282	chuckv	1388	real(kind = dp), intent(inout) :: rho_i_at_j
283	chuckv	702	! value of electron density rho do to atom j at atom i
284	chuckv	1388	real(kind = dp), intent(inout) :: rho_j_at_i
285	chuckv	702	! we don't use the derivatives, dummy variables
286	chuckv	728	real( kind = dp) :: drho
287	chuckv	714	integer :: sc_err
288	chuckv	702
289	gezelter	1285	integer :: myid_atom1 ! SC atid
290	chuckv	702	integer :: myid_atom2
291	gezelter	1448
292	chuckv	702	! check to see if we need to be cleaned at the start of a force loop
293	gezelter	1448
294	chuckv	1388	if (.not.haveMixingMap) call createMixingMap()
295			haveMixingMap = .true.
296	chuckv	1175
297	chuckv	728	Myid_atom1 = SCList%atidtoSCtype(Atid1)
298			Myid_atom2 = SCList%atidtoSCtype(Atid2)
299	gezelter	938
300			call lookupUniformSpline(MixingMap(myid_atom1,myid_atom2)%phi, r, &
301			rho_i_at_j)
302			rho_j_at_i = rho_i_at_j
303	chuckv	1388
304	chuckv	713	end subroutine calc_sc_prepair_rho
305	gezelter	1479
306	gezelter	938	!! Calculate the rho_a for all local atoms
307	gezelter	1479	subroutine calc_sc_preforce_Frho(atid, rho, frho, dfrhodrho)
308			real(kind=dp) :: rho, frho, dfrhodrho
309			real(kind=dp) :: u, u1
310			integer :: atid, myid
311	chuckv	1175
312	chuckv	713	!! Calculate F(rho) and derivative
313	gezelter	1479
314			Myid = SCList%atidtoSctype(atid)
315			! Myid is set to -1 for non SC atoms.
316			! Punt if we are a non-SC atom type.
317			if (Myid == -1) then
318			frho = 0.0_dp
319			dfrhodrho = 0.0_dp
320			else
321			frho = - SCList%SCTypes(Myid)%c * &
322			SCList%SCTypes(Myid)%epsilon * sqrt(rho)
323	gezelter	1448
324	gezelter	1479	dfrhodrho = 0.5_dp*frho/rho
325			end if
326	gezelter	1448
327			end subroutine calc_sc_preforce_Frho
328	gezelter	938
329	chuckv	713	!! Does Sutton-Chen pairwise Force calculation.
330	gezelter	1464	subroutine do_sc_pair(atid1, atid2, d, rij, r2, sw, vpair, &
331	gezelter	1467	pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, &
332	gezelter	1464	fshift_i, fshift_j)
333	chuckv	702	!Arguments
334	gezelter	1464	integer, intent(in) :: atid1, atid2
335	gezelter	1419	real( kind = dp ), intent(in) :: rij, r2
336	chuckv	702	real( kind = dp ) :: pot, sw, vpair
337	gezelter	1386	real( kind = dp ), dimension(3) :: f1
338	chuckv	702	real( kind = dp ), intent(in), dimension(3) :: d
339	chuckv	1388	real( kind = dp ), intent(inout) :: dfrhodrho_i, dfrhodrho_j
340			real( kind = dp ), intent(inout) :: rho_i, rho_j
341	gezelter	1464	real( kind = dp ), intent(inout):: fshift_i, fshift_j
342	gezelter	1448
343	gezelter	938	real( kind = dp ) :: drdx, drdy, drdz
344	chuckv	728	real( kind = dp ) :: dvpdr
345	gezelter	938	real( kind = dp ) :: rhtmp, drhodr
346	chuckv	702	real( kind = dp ) :: dudr
347			real( kind = dp ) :: Fx,Fy,Fz
348	gezelter	938	real( kind = dp ) :: pot_temp, vptmp
349			real( kind = dp ) :: rcij, vcij
350	gezelter	1448
351	gezelter	938	integer :: id1, id2
352	chuckv	702	integer :: mytype_atom1
353			integer :: mytype_atom2
354	gezelter	1448
355	gezelter	1386	mytype_atom1 = SCList%atidToSCType(atid1)
356			mytype_atom2 = SCList%atidTOSCType(atid2)
357
358			drdx = d(1)/rij
359			drdy = d(2)/rij
360			drdz = d(3)/rij
361
362			rcij = MixingMap(mytype_atom1,mytype_atom2)%rCut
363			vcij = MixingMap(mytype_atom1,mytype_atom2)%vCut
364
365			call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%phi, &
366			rij, rhtmp, drhodr)
367
368			call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%V, &
369			rij, vptmp, dvpdr)
370	gezelter	1448
371	chuckv	1388	dudr = drhodr*(dfrhodrho_i + dfrhodrho_j) + dvpdr
372
373			pot_temp = vptmp - vcij
374	gezelter	1386
375	chuckv	1388	vpair = vpair + pot_temp
376	gezelter	1448
377	chuckv	1388	fx = dudr * drdx
378			fy = dudr * drdy
379			fz = dudr * drdz
380	gezelter	1464
381			! particle_pot is the difference between the full potential
382			! and the full potential without the presence of a particular
383			! particle (atom1).
384			!
385			! This reduces the density at other particle locations, so
386			! we need to recompute the density at atom2 assuming atom1
387			! didn't contribute. This then requires recomputing the
388			! density functional for atom2 as well.
389			!
390			! Most of the particle_pot heavy lifting comes from the
391			! pair interaction, and will be handled by vpair.
392	gezelter	1448
393	gezelter	1464	fshift_i = -SCList%SCTypes(mytype_atom1)%c * &
394			SCList%SCTypes(mytype_atom1)%epsilon * sqrt(rho_i-rhtmp)
395			fshift_j = -SCList%SCTypes(mytype_atom2)%c * &
396			SCList%SCTypes(mytype_atom2)%epsilon * sqrt(rho_j-rhtmp)
397	gezelter	1448
398	chuckv	1388	pot = pot + pot_temp
399
400			f1(1) = f1(1) + fx
401			f1(2) = f1(2) + fy
402			f1(3) = f1(3) + fz
403	gezelter	1448
404	chuckv	728	end subroutine do_sc_pair
405	chuckv	713	end module suttonchen