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(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(atid1, atid2, d, rij, r2, sw, vpair, &
       fpair, 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), 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   
    
    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:	1465
Committed:	Fri Jul 9 23:08:25 2010 UTC (15 years ago) by chuckv
File size:	13471 byte(s)
Log Message:	Creating busticated version of OpenMD
#	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	1448
306
307	gezelter	938	!! Calculate the rho_a for all local atoms
308	gezelter	1285	subroutine calc_sc_preforce_Frho(nlocal, pot, particle_pot)
309	chuckv	702	integer :: nlocal
310			real(kind=dp) :: pot
311	gezelter	1285	real(kind=dp), dimension(nlocal) :: particle_pot
312	chuckv	702	integer :: i,j
313			integer :: atom
314			integer :: atype1
315	chuckv	728	integer :: atid1
316			integer :: myid
317	chuckv	1175
318	chuckv	713	!! Calculate F(rho) and derivative
319	chuckv	702	do atom = 1, nlocal
320	chuckv	728	Myid = SCList%atidtoSctype(Atid(atom))
321	chuckv	1175	! Myid is set to -1 for non SC atoms.
322			! Punt if we are a non-SC atom type.
323			if (Myid == -1) then
324			frho(atom) = 0.0_dp
325			dfrhodrho(atom) = 0.0_dp
326			else
327			frho(atom) = - SCList%SCTypes(Myid)%c * &
328	gezelter	1448	SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom))
329
330	chuckv	1175	dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
331			end if
332	chuckv	842	pot = pot + frho(atom)
333	gezelter	1285	particle_pot(atom) = particle_pot(atom) + frho(atom)
334	chuckv	702	enddo
335	gezelter	1448
336			end subroutine calc_sc_preforce_Frho
337	gezelter	938
338	chuckv	713	!! Does Sutton-Chen pairwise Force calculation.
339	gezelter	1464	subroutine do_sc_pair(atid1, atid2, d, rij, r2, sw, vpair, &
340	gezelter	1448	fpair, pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, &
341	gezelter	1464	fshift_i, fshift_j)
342	chuckv	702	!Arguments
343	gezelter	1464	integer, intent(in) :: atid1, atid2
344	gezelter	1419	real( kind = dp ), intent(in) :: rij, r2
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	gezelter	1464	real( kind = dp ), intent(inout):: fshift_i, fshift_j
352	gezelter	1448
353	gezelter	938	real( kind = dp ) :: drdx, drdy, drdz
354	chuckv	728	real( kind = dp ) :: dvpdr
355	gezelter	938	real( kind = dp ) :: rhtmp, drhodr
356	chuckv	702	real( kind = dp ) :: dudr
357			real( kind = dp ) :: Fx,Fy,Fz
358	gezelter	938	real( kind = dp ) :: pot_temp, vptmp
359			real( kind = dp ) :: rcij, vcij
360	gezelter	1448
361	gezelter	938	integer :: id1, id2
362	chuckv	702	integer :: mytype_atom1
363			integer :: mytype_atom2
364	gezelter	1448
365	gezelter	1386	mytype_atom1 = SCList%atidToSCType(atid1)
366			mytype_atom2 = SCList%atidTOSCType(atid2)
367
368			drdx = d(1)/rij
369			drdy = d(2)/rij
370			drdz = d(3)/rij
371
372			rcij = MixingMap(mytype_atom1,mytype_atom2)%rCut
373			vcij = MixingMap(mytype_atom1,mytype_atom2)%vCut
374
375			call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%phi, &
376			rij, rhtmp, drhodr)
377
378			call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%V, &
379			rij, vptmp, dvpdr)
380	gezelter	1448
381	chuckv	1388	dudr = drhodr*(dfrhodrho_i + dfrhodrho_j) + dvpdr
382
383			pot_temp = vptmp - vcij
384	gezelter	1386
385	chuckv	1388	vpair = vpair + pot_temp
386	gezelter	1448
387	chuckv	1388	fx = dudr * drdx
388			fy = dudr * drdy
389			fz = dudr * drdz
390	gezelter	1464
391			! particle_pot is the difference between the full potential
392			! and the full potential without the presence of a particular
393			! particle (atom1).
394			!
395			! This reduces the density at other particle locations, so
396			! we need to recompute the density at atom2 assuming atom1
397			! didn't contribute. This then requires recomputing the
398			! density functional for atom2 as well.
399			!
400			! Most of the particle_pot heavy lifting comes from the
401			! pair interaction, and will be handled by vpair.
402	gezelter	1448
403	gezelter	1464	fshift_i = -SCList%SCTypes(mytype_atom1)%c * &
404			SCList%SCTypes(mytype_atom1)%epsilon * sqrt(rho_i-rhtmp)
405			fshift_j = -SCList%SCTypes(mytype_atom2)%c * &
406			SCList%SCTypes(mytype_atom2)%epsilon * sqrt(rho_j-rhtmp)
407	gezelter	1448
408	chuckv	1388	pot = pot + pot_temp
409
410			f1(1) = f1(1) + fx
411			f1(2) = f1(2) + fy
412			f1(3) = f1(3) + fz
413	gezelter	1448
414	chuckv	728	end subroutine do_sc_pair
415	chuckv	713	end module suttonchen