ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/eam.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/eam.F90 (file contents):
Revision 1624 by chuckv, Thu Oct 21 15:25:30 2004 UTC vs.
Revision 2728 by chrisfen, Fri Apr 21 20:02:19 2006 UTC

# Line 1 | Line 1
1 + !!
2 + !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 + !!
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   module eam
43    use simulation
44    use force_globals
45    use status
46    use atype_module
47 <  use Vector_class
47 >  use vector_class
48 >  use interpolation
49   #ifdef IS_MPI
50    use mpiSimulation
51   #endif
52    implicit none
53    PRIVATE
54 + #define __FORTRAN90
55 + #include "UseTheForce/DarkSide/fInteractionMap.h"
56  
57    INTEGER, PARAMETER :: DP = selected_real_kind(15)
58  
# Line 22 | Line 66 | module eam
66  
67    character(len = 200) :: errMsg
68    character(len=*), parameter :: RoutineName =  "EAM MODULE"
69 < !! Logical that determines if eam arrays should be zeroed
69 >  !! Logical that determines if eam arrays should be zeroed
70    logical :: cleanme = .true.
27 –  logical :: nmflag  = .false.
71  
29 –
72    type, private :: EAMtype
73       integer           :: eam_atype      
32 –     real( kind = DP ) :: eam_dr          
33 –     integer           :: eam_nr          
34 –     integer           :: eam_nrho          
74       real( kind = DP ) :: eam_lattice        
36 –     real( kind = DP ) :: eam_drho      
75       real( kind = DP ) :: eam_rcut    
76       integer           :: eam_atype_map
77 <    
78 <     real( kind = DP ), pointer, dimension(:) :: eam_rvals        => null()
79 <     real( kind = DP ), pointer, dimension(:) :: eam_rhovals      => null()
80 <     real( kind = DP ), pointer, dimension(:) :: eam_F_rho        => null()
43 <     real( kind = DP ), pointer, dimension(:) :: eam_Z_r          => null()
44 <     real( kind = DP ), pointer, dimension(:) :: eam_rho_r        => null()
45 <     real( kind = DP ), pointer, dimension(:) :: eam_phi_r        => null()
46 <     real( kind = DP ), pointer, dimension(:) :: eam_F_rho_pp     => null()
47 <     real( kind = DP ), pointer, dimension(:) :: eam_Z_r_pp       => null()
48 <     real( kind = DP ), pointer, dimension(:) :: eam_rho_r_pp     => null()
49 <     real( kind = DP ), pointer, dimension(:) :: eam_phi_r_pp     => null()
77 >     type(cubicSpline) :: rho
78 >     type(cubicSpline) :: Z
79 >     type(cubicSpline) :: F
80 >     type(cubicSpline) :: phi
81    end type EAMtype
82  
52 –
83    !! Arrays for derivatives used in force calculation
84    real( kind = dp), dimension(:), allocatable :: frho
85    real( kind = dp), dimension(:), allocatable :: rho
56 –
86    real( kind = dp), dimension(:), allocatable :: dfrhodrho
58 –  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
87  
88 <
61 < !! Arrays for MPI storage
88 >  !! Arrays for MPI storage
89   #ifdef IS_MPI
90    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
91    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
# Line 67 | Line 94 | module eam
94    real( kind = dp),save, dimension(:), allocatable :: rho_row
95    real( kind = dp),save, dimension(:), allocatable :: rho_col
96    real( kind = dp),save, dimension(:), allocatable :: rho_tmp
70 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
71 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
97   #endif
98  
99    type, private :: EAMTypeList
100       integer           :: n_eam_types = 0
101       integer           :: currentAddition = 0
102 <    
102 >
103       type (EAMtype), pointer  :: EAMParams(:) => null()
104 +     integer, pointer         :: atidToEAMType(:) => null()
105    end type EAMTypeList
106  
81 –
107    type (eamTypeList), save :: EAMList
108  
109    !! standard eam stuff  
# Line 91 | Line 116 | module eam
116    public :: calc_eam_prepair_rho
117    public :: calc_eam_preforce_Frho
118    public :: clean_EAM
119 +  public :: destroyEAMTypes
120 +  public :: getEAMCut
121 +  public :: lookupEAMSpline
122 +  public :: lookupEAMSpline1d
123  
124   contains
125  
97 –
126    subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
127 <       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
100 <       eam_ident,status)
127 >       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho, c_ident, status)
128      real (kind = dp )                      :: lattice_constant
129      integer                                :: eam_nrho
130      real (kind = dp )                      :: eam_drho
131      integer                                :: eam_nr
132      real (kind = dp )                      :: eam_dr
133      real (kind = dp )                      :: rcut
134 <    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r
135 <    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
136 <    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
137 <    integer                                :: eam_ident
134 >    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r, rvals
135 >    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r, eam_phi_r
136 >    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho, rhovals
137 >    integer                                :: c_ident
138      integer                                :: status
139  
140 <    integer                                :: nAtypes
140 >    integer                                :: nAtypes,nEAMTypes,myATID
141      integer                                :: maxVals
142      integer                                :: alloc_stat
143 <    integer                                :: current
143 >    integer                                :: current, j
144      integer,pointer                        :: Matchlist(:) => null()
145  
146      status = 0
147  
121 –
148      !! Assume that atypes has already been set and get the total number of types in atypes
149      !! Also assume that every member of atypes is a EAM model.
124 –  
150  
151      ! check to see if this is the first time into
152      if (.not.associated(EAMList%EAMParams)) then
153 <       call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
154 <       EAMList%n_eam_types = nAtypes
155 <       allocate(EAMList%EAMParams(nAtypes))
153 >       call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
154 >       EAMList%n_eam_types = nEAMtypes
155 >       allocate(EAMList%EAMParams(nEAMTypes))
156 >       nAtypes = getSize(atypes)
157 >       allocate(EAMList%atidToEAMType(nAtypes))
158      end if
159  
160      EAMList%currentAddition = EAMList%currentAddition + 1
161      current = EAMList%currentAddition
135 –    
162  
163 <    call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
164 <    if (alloc_stat /= 0) then
139 <       status = -1
140 <       return
141 <    end if
163 >    myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
164 >    EAMList%atidToEAMType(myATID) = current
165  
166 <    ! this is a possible bug, we assume a correspondence between the vector atypes and
144 <    ! EAMAtypes
145 <      
146 <    EAMList%EAMParams(current)%eam_atype    = eam_ident
166 >    EAMList%EAMParams(current)%eam_atype    = c_ident
167      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
148 –    EAMList%EAMParams(current)%eam_nrho     = eam_nrho
149 –    EAMList%EAMParams(current)%eam_drho     = eam_drho
150 –    EAMList%EAMParams(current)%eam_nr       = eam_nr
151 –    EAMList%EAMParams(current)%eam_dr       = eam_dr
168      EAMList%EAMParams(current)%eam_rcut     = rcut
153 –    EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
154 –    EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
155 –    EAMList%EAMParams(current)%eam_F_rho    = eam_F_rho
169  
170 +    ! Build array of r values
171 +    do j = 1, eam_nr
172 +       rvals(j) = real(j-1,kind=dp) * eam_dr
173 +    end do
174 +    
175 +    ! Build array of rho values
176 +    do j = 1, eam_nrho
177 +       rhovals(j) = real(j-1,kind=dp) * eam_drho
178 +    end do
179 +
180 +    ! convert from eV to kcal / mol:
181 +    do j = 1, eam_nrho
182 +       eam_F_rho(j) = eam_F_rho(j) * 23.06054E0_DP
183 +    end do
184 +    
185 +    ! precompute the pair potential and get it into kcal / mol:
186 +    eam_phi_r(1) = 0.0E0_DP
187 +    do j = 2, eam_nr
188 +       eam_phi_r(j) = 331.999296E0_DP * (eam_Z_r(j)**2) / rvals(j)
189 +    enddo
190 +
191 +    call newSpline(EAMList%EAMParams(current)%rho, rvals, eam_rho_r, .true.)
192 +    call newSpline(EAMList%EAMParams(current)%Z,   rvals, eam_Z_r, .true.)
193 +    call newSpline(EAMList%EAMParams(current)%F, rhovals, eam_F_rho, .true.)
194 +    call newSpline(EAMList%EAMParams(current)%phi, rvals, eam_phi_r, .true.)
195    end subroutine newEAMtype
196  
197  
198 +  ! kills all eam types entered and sets simulation to uninitalized
199 +  subroutine destroyEAMtypes()
200 +    integer :: i
201 +    type(EAMType), pointer :: tempEAMType=>null()
202  
203 +    do i = 1, EAMList%n_eam_types
204 +       tempEAMType => eamList%EAMParams(i)
205 +       call deallocate_EAMType(tempEAMType)
206 +    end do
207 +    if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
208 +    eamList%EAMParams => null()
209 +
210 +    eamList%n_eam_types = 0
211 +    eamList%currentAddition = 0
212 +  end subroutine destroyEAMtypes
213 +
214 +  function getEAMCut(atomID) result(cutValue)
215 +    integer, intent(in) :: atomID
216 +    integer :: eamID
217 +    real(kind=dp) :: cutValue
218 +    
219 +    eamID = EAMList%atidToEAMType(atomID)
220 +    cutValue = EAMList%EAMParams(eamID)%eam_rcut
221 +  end function getEAMCut
222 +
223    subroutine init_EAM_FF(status)
224      integer :: status
225      integer :: i,j
226      real(kind=dp) :: current_rcut_max
227 + #ifdef IS_MPI
228 +    integer :: nAtomsInRow
229 +    integer :: nAtomsInCol
230 + #endif
231      integer :: alloc_stat
232      integer :: number_r, number_rho
233  
168 –
234      status = 0
235      if (EAMList%currentAddition == 0) then
236         call handleError("init_EAM_FF","No members in EAMList")
237         status = -1
238         return
239      end if
240 <
241 <
177 <       do i = 1, EAMList%currentAddition
178 <
179 < ! Build array of r values
180 <
181 <          do j = 1,EAMList%EAMParams(i)%eam_nr
182 <             EAMList%EAMParams(i)%eam_rvals(j) = &
183 <                  real(j-1,kind=dp)* &
184 <                  EAMList%EAMParams(i)%eam_dr
185 <              end do
186 < ! Build array of rho values
187 <          do j = 1,EAMList%EAMParams(i)%eam_nrho
188 <             EAMList%EAMParams(i)%eam_rhovals(j) = &
189 <                  real(j-1,kind=dp)* &
190 <                  EAMList%EAMParams(i)%eam_drho
191 <          end do
192 <          ! convert from eV to kcal / mol:
193 <          EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
194 <
195 <          ! precompute the pair potential and get it into kcal / mol:
196 <          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
197 <          do j = 2, EAMList%EAMParams(i)%eam_nr
198 <             EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
199 <             EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
200 <          enddo
201 <       end do
202 <      
203 <
204 <       do i = 1,  EAMList%currentAddition
205 <          number_r   = EAMList%EAMParams(i)%eam_nr
206 <          number_rho = EAMList%EAMParams(i)%eam_nrho
207 <          
208 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
209 <               EAMList%EAMParams(i)%eam_rho_r, &
210 <               EAMList%EAMParams(i)%eam_rho_r_pp, &
211 <               0.0E0_DP, 0.0E0_DP, 'N')
212 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
213 <               EAMList%EAMParams(i)%eam_Z_r, &
214 <               EAMList%EAMParams(i)%eam_Z_r_pp, &
215 <               0.0E0_DP, 0.0E0_DP, 'N')
216 <          call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
217 <               EAMList%EAMParams(i)%eam_F_rho, &
218 <               EAMList%EAMParams(i)%eam_F_rho_pp, &
219 <               0.0E0_DP, 0.0E0_DP, 'N')
220 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
221 <               EAMList%EAMParams(i)%eam_phi_r, &
222 <               EAMList%EAMParams(i)%eam_phi_r_pp, &
223 <               0.0E0_DP, 0.0E0_DP, 'N')
224 <       enddo
225 <
226 < !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
227 <       !! find the smallest rcut for any eam atype
228 < !       do i = 2, EAMList%currentAddition
229 < !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
230 < !       end do
231 <
232 < !       EAM_rcut = current_rcut_max
233 < !       EAM_rcut_orig = current_rcut_max
234 < !       do i = 1, EAMList%currentAddition
235 < !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
236 < !       end do
237 <       !! Allocate arrays for force calculation
238 <      
239 <       call allocateEAM(alloc_stat)
240 <       if (alloc_stat /= 0 ) then
241 <          write(*,*) "allocateEAM failed"
242 <          status = -1
243 <          return
244 <       endif
240 >    
241 >    !! Allocate arrays for force calculation
242      
246 –  end subroutine init_EAM_FF
247 –
248 – !! routine checks to see if array is allocated, deallocates array if allocated
249 – !! and then creates the array to the required size
250 –  subroutine allocateEAM(status)
251 –    integer, intent(out) :: status
252 –
243   #ifdef IS_MPI
254 –    integer :: nAtomsInRow
255 –    integer :: nAtomsInCol
256 – #endif
257 –    integer :: alloc_stat
258 –
259 –
260 –    status = 0
261 – #ifdef IS_MPI
244      nAtomsInRow = getNatomsInRow(plan_atom_row)
245      nAtomsInCol = getNatomsInCol(plan_atom_col)
246   #endif
# Line 269 | Line 251 | contains
251         status = -1
252         return
253      end if
254 +
255      if (allocated(rho)) deallocate(rho)
256      allocate(rho(nlocal),stat=alloc_stat)
257      if (alloc_stat /= 0) then
# Line 283 | Line 266 | contains
266         return
267      end if
268  
286 –    if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
287 –    allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
288 –    if (alloc_stat /= 0) then
289 –       status = -1
290 –       return
291 –    end if
292 –    
269   #ifdef IS_MPI
270  
271      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 299 | Line 275 | contains
275         return
276      end if
277  
302 –
278      if (allocated(frho_row)) deallocate(frho_row)
279      allocate(frho_row(nAtomsInRow),stat=alloc_stat)
280      if (alloc_stat /= 0) then
# Line 318 | Line 293 | contains
293         status = -1
294         return
295      end if
321 –    if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
322 –    allocate(d2frhodrhodrho_row(nAtomsInRow),stat=alloc_stat)
323 –    if (alloc_stat /= 0) then
324 –       status = -1
325 –       return
326 –    end if
296  
297 +    ! Now do column arrays
298  
329 – ! Now do column arrays
330 –
299      if (allocated(frho_col)) deallocate(frho_col)
300      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
301      if (alloc_stat /= 0) then
# Line 346 | Line 314 | contains
314         status = -1
315         return
316      end if
317 <    if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
350 <    allocate(d2frhodrhodrho_col(nAtomsInCol),stat=alloc_stat)
351 <    if (alloc_stat /= 0) then
352 <       status = -1
353 <       return
354 <    end if
355 <  
317 >
318   #endif
319  
320 <  end subroutine allocateEAM
320 >  end subroutine init_EAM_FF
321  
322 < !! C sets rcut to be the largest cutoff of any atype
361 < !! present in this simulation. Doesn't include all atypes
362 < !! sim knows about, just those in the simulation.
363 <  subroutine setCutoffEAM(rcut, status)
322 >  subroutine setCutoffEAM(rcut)
323      real(kind=dp) :: rcut
365 –    integer :: status
366 –    status = 0
367 –
324      EAM_rcut = rcut
369 –
325    end subroutine setCutoffEAM
326  
372 –
373 –
327    subroutine clean_EAM()
328 <  
329 <   ! clean non-IS_MPI first
328 >
329 >    ! clean non-IS_MPI first
330      frho = 0.0_dp
331      rho  = 0.0_dp
332      dfrhodrho = 0.0_dp
333 < ! clean MPI if needed
333 >    ! clean MPI if needed
334   #ifdef IS_MPI
335      frho_row = 0.0_dp
336      frho_col = 0.0_dp
# Line 389 | Line 342 | contains
342   #endif
343    end subroutine clean_EAM
344  
392 –
393 –
394 –  subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
395 –    integer, intent(in)          :: eam_n_rho
396 –    integer, intent(in)          :: eam_n_r
397 –    type (EAMType)               :: thisEAMType
398 –    integer, optional   :: stat
399 –    integer             :: alloc_stat
400 –
401 –
402 –
403 –    if (present(stat)) stat = 0
404 –    
405 –    allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
406 –    if (alloc_stat /= 0 ) then
407 –       if (present(stat)) stat = -1
408 –       return
409 –    end if
410 –    allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat)  
411 –    if (alloc_stat /= 0 ) then
412 –       if (present(stat)) stat = -1
413 –       return
414 –    end if
415 –    allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)  
416 –    if (alloc_stat /= 0 ) then
417 –       if (present(stat)) stat = -1
418 –       return
419 –    end if
420 –    allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)        
421 –    if (alloc_stat /= 0 ) then
422 –       if (present(stat)) stat = -1
423 –       return
424 –    end if
425 –    allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat)      
426 –    if (alloc_stat /= 0 ) then
427 –       if (present(stat)) stat = -1
428 –       return
429 –    end if
430 –    allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat)      
431 –    if (alloc_stat /= 0 ) then
432 –       if (present(stat)) stat = -1
433 –       return
434 –    end if
435 –    allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat)  
436 –    if (alloc_stat /= 0 ) then
437 –       if (present(stat)) stat = -1
438 –       return
439 –    end if
440 –    allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)  
441 –    if (alloc_stat /= 0 ) then
442 –       if (present(stat)) stat = -1
443 –       return
444 –    end if
445 –    allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)  
446 –    if (alloc_stat /= 0 ) then
447 –       if (present(stat)) stat = -1
448 –       return
449 –    end if
450 –    allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
451 –    if (alloc_stat /= 0 ) then
452 –       if (present(stat)) stat = -1
453 –       return
454 –    end if
455 –      
456 –
457 –  end subroutine allocate_EAMType
458 –
459 –
345    subroutine deallocate_EAMType(thisEAMType)
346      type (EAMtype), pointer :: thisEAMType
347  
348 <    ! free Arrays in reverse order of allocation...
349 <    deallocate(thisEAMType%eam_phi_r_pp)      
350 <    deallocate(thisEAMType%eam_rho_r_pp)  
351 <    deallocate(thisEAMType%eam_Z_r_pp)  
352 <    deallocate(thisEAMType%eam_F_rho_pp)  
468 <    deallocate(thisEAMType%eam_phi_r)      
469 <    deallocate(thisEAMType%eam_rho_r)      
470 <    deallocate(thisEAMType%eam_Z_r)  
471 <    deallocate(thisEAMType%eam_F_rho)
472 <    deallocate(thisEAMType%eam_rhovals)
473 <    deallocate(thisEAMType%eam_rvals)
474 <  
348 >    call deleteSpline(thisEAMType%F)
349 >    call deleteSpline(thisEAMType%rho)
350 >    call deleteSpline(thisEAMType%phi)
351 >    call deleteSpline(thisEAMType%Z)
352 >
353    end subroutine deallocate_EAMType
354  
355 < !! Calculates rho_r
355 >  !! Calculates rho_r
356    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
357 <    integer :: atom1,atom2
357 >    integer :: atom1, atom2
358      real(kind = dp), dimension(3) :: d
359      real(kind = dp), intent(inout)               :: r
360      real(kind = dp), intent(inout)               :: rijsq
# Line 484 | Line 362 | contains
362      real(kind = dp) :: rho_i_at_j
363      ! value of electron density rho do to atom j at atom i
364      real(kind = dp) :: rho_j_at_i
487 –
488 –    ! we don't use the derivatives, dummy variables
489 –    real( kind = dp) :: drho,d2rho
365      integer :: eam_err
366 <  
367 <    integer :: myid_atom1
368 <    integer :: myid_atom2
366 >    
367 >    integer :: atid1, atid2 ! Global atid    
368 >    integer :: myid_atom1 ! EAM atid
369 >    integer :: myid_atom2
370  
371 < ! check to see if we need to be cleaned at the start of a force loop
496 <      
371 >    ! check to see if we need to be cleaned at the start of a force loop
372  
498 –
499 –
373   #ifdef IS_MPI
374 <    myid_atom1 = atid_Row(atom1)
375 <    myid_atom2 = atid_Col(atom2)
374 >    Atid1 = Atid_row(Atom1)
375 >    Atid2 = Atid_col(Atom2)
376   #else
377 <    myid_atom1 = atid(atom1)
378 <    myid_atom2 = atid(atom2)
377 >    Atid1 = Atid(Atom1)
378 >    Atid2 = Atid(Atom2)
379   #endif
380  
381 +    Myid_atom1 = Eamlist%atidtoeamtype(Atid1)
382 +    Myid_atom2 = Eamlist%atidtoeamtype(Atid2)
383 +
384      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
385  
386 +       call lookupEAMSpline(EAMList%EAMParams(myid_atom1)%rho, r, &
387 +            rho_i_at_j)
388  
511 –
512 –       call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
513 –            EAMList%EAMParams(myid_atom1)%eam_rvals, &
514 –            EAMList%EAMParams(myid_atom1)%eam_rho_r, &
515 –            EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
516 –            r, rho_i_at_j,drho,d2rho)
517 –
518 –
519 –      
389   #ifdef  IS_MPI
390         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
391   #else
392         rho(atom2) = rho(atom2) + rho_i_at_j
393   #endif
394 < !       write(*,*) atom1,atom2,r,rho_i_at_j
526 <       endif
394 >    endif
395  
396 <       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
529 <          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
530 <               EAMList%EAMParams(myid_atom2)%eam_rvals, &
531 <               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
532 <               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
533 <               r, rho_j_at_i,drho,d2rho)
396 >    if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
397  
398 <    
399 <      
400 <      
398 >       call lookupEAMSpline(EAMList%EAMParams(myid_atom2)%rho, r, &
399 >            rho_j_at_i)
400 >
401   #ifdef  IS_MPI
402 <          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
402 >       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
403   #else
404 <          rho(atom1) = rho(atom1) + rho_j_at_i
404 >       rho(atom1) = rho(atom1) + rho_j_at_i
405   #endif
406 <       endif
544 <
545 <
546 <
547 <
548 <
549 <
406 >    endif
407    end subroutine calc_eam_prepair_rho
408  
409  
553 –
554 –
410    !! Calculate the functional F(rho) for all local atoms
411 <  subroutine calc_eam_preforce_Frho(nlocal,pot)
411 >  subroutine calc_eam_preforce_Frho(nlocal, pot)
412      integer :: nlocal
413      real(kind=dp) :: pot
414 <    integer :: i,j
414 >    integer :: i, j
415      integer :: atom
416 <    real(kind=dp) :: U,U1,U2
416 >    real(kind=dp) :: U,U1
417      integer :: atype1
418 <    integer :: me
564 <    integer :: n_rho_points
418 >    integer :: me, atid1
419  
566 –  
420      cleanme = .true.
421 <    !! Scatter the electron density from  pre-pair calculation back to local atoms
421 >    !! Scatter the electron density from  pre-pair calculation back to
422 >    !! local atoms
423   #ifdef IS_MPI
424      call scatter(rho_row,rho,plan_atom_row,eam_err)
425      if (eam_err /= 0 ) then
426 <      write(errMsg,*) " Error scattering rho_row into rho"
427 <      call handleError(RoutineName,errMesg)
428 <   endif      
426 >       write(errMsg,*) " Error scattering rho_row into rho"
427 >       call handleError(RoutineName,errMesg)
428 >    endif
429      call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
430      if (eam_err /= 0 ) then
431 <      write(errMsg,*) " Error scattering rho_col into rho"
432 <      call handleError(RoutineName,errMesg)
433 <   endif
431 >       write(errMsg,*) " Error scattering rho_col into rho"
432 >       call handleError(RoutineName,errMesg)
433 >    endif
434  
435 <      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
435 >    rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
436   #endif
437  
438 <
585 <
586 < !! Calculate F(rho) and derivative
438 >    !! Calculate F(rho) and derivative
439      do atom = 1, nlocal
440 <       me = atid(atom)
441 <       n_rho_points = EAMList%EAMParams(me)%eam_nrho
590 <       !  Check to see that the density is not greater than the larges rho we have calculated
591 <       if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
592 <          call eam_splint(n_rho_points, &
593 <               EAMList%EAMParams(me)%eam_rhovals, &
594 <               EAMList%EAMParams(me)%eam_f_rho, &
595 <               EAMList%EAMParams(me)%eam_f_rho_pp, &
596 <               rho(atom), & ! Actual Rho
597 <               u, u1, u2)
598 <       else
599 <          ! Calculate F(rho with the largest available rho value
600 <          call eam_splint(n_rho_points, &
601 <               EAMList%EAMParams(me)%eam_rhovals, &
602 <               EAMList%EAMParams(me)%eam_f_rho, &
603 <               EAMList%EAMParams(me)%eam_f_rho_pp, &
604 <               EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
605 <               u,u1,u2)
606 <       end if
440 >       atid1 = atid(atom)
441 >       me = eamList%atidToEAMtype(atid1)
442  
443 <
443 >       call lookupEAMSpline1d(EAMList%EAMParams(me)%F, rho(atom), &
444 >            u, u1)
445 >      
446         frho(atom) = u
447         dfrhodrho(atom) = u1
611 –       d2frhodrhodrho(atom) = u2
448         pot = pot + u
449  
450      enddo
451  
616 –  
617 –
452   #ifdef IS_MPI
453      !! communicate f(rho) and derivatives back into row and column arrays
454      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 633 | Line 467 | contains
467      if (eam_err /=  0) then
468         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
469      endif
636 –
637 –
638 –
639 –
640 –
641 –    if (nmflag) then
642 –       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row)
643 –       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col)
644 –    endif
470   #endif
471  
647 –  
648 –  end subroutine calc_eam_preforce_Frho
472  
473 <
474 <
652 <
473 >  end subroutine calc_eam_preforce_Frho
474 >  
475    !! Does EAM pairwise Force calculation.  
476    subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
477         pot, f, do_pot)
# Line 662 | Line 484 | contains
484      real( kind = dp ), intent(inout), dimension(3) :: fpair
485  
486      logical, intent(in) :: do_pot
487 <    
488 <    real( kind = dp ) :: drdx,drdy,drdz
489 <    real( kind = dp ) :: d2
490 <    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
491 <    real( kind = dp ) :: rha,drha,d2rha, dpha
670 <    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
487 >
488 >    real( kind = dp ) :: drdx, drdy, drdz
489 >    real( kind = dp ) :: phab, pha, dvpdr
490 >    real( kind = dp ) :: rha, drha, dpha
491 >    real( kind = dp ) :: rhb, drhb, dphb
492      real( kind = dp ) :: dudr
493 <    real( kind = dp ) :: rci,rcj
494 <    real( kind = dp ) :: drhoidr,drhojdr
495 <    real( kind = dp ) :: d2rhoidrdr
496 <    real( kind = dp ) :: d2rhojdrdr
676 <    real( kind = dp ) :: Fx,Fy,Fz
677 <    real( kind = dp ) :: r,d2pha,phb,d2phb
493 >    real( kind = dp ) :: rci, rcj
494 >    real( kind = dp ) :: drhoidr, drhojdr
495 >    real( kind = dp ) :: Fx, Fy, Fz
496 >    real( kind = dp ) :: r, phb
497  
498 <    integer :: id1,id2
498 >    integer :: id1, id2
499      integer  :: mytype_atom1
500      integer  :: mytype_atom2
501 +    integer  :: atid1, atid2
502  
683 – !Local Variables
684 –    
685 –   ! write(*,*) "Frho: ", Frho(atom1)
686 –
503      phab = 0.0E0_DP
504      dvpdr = 0.0E0_DP
689 –    d2vpdrdr = 0.0E0_DP
505  
506      if (rij .lt. EAM_rcut) then
507  
508   #ifdef IS_MPI
509 <       mytype_atom1 = atid_row(atom1)
510 <       mytype_atom2 = atid_col(atom2)
509 >       atid1 = atid_row(atom1)
510 >       atid2 = atid_col(atom2)
511   #else
512 <       mytype_atom1 = atid(atom1)
513 <       mytype_atom2 = atid(atom2)
512 >       atid1 = atid(atom1)
513 >       atid2 = atid(atom2)
514   #endif
515 +
516 +       mytype_atom1 = EAMList%atidToEAMType(atid1)
517 +       mytype_atom2 = EAMList%atidTOEAMType(atid2)
518 +
519 +
520         ! get cutoff for atom 1
521         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
522         ! get type specific cutoff for atom 2
523         rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
524 <      
524 >
525         drdx = d(1)/rij
526         drdy = d(2)/rij
527         drdz = d(3)/rij
528 <      
528 >
529         if (rij.lt.rci) then
530 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
531 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
532 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
533 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
534 <               rij, rha,drha,d2rha)
535 <          !! Calculate Phi(r) for atom1.
536 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
537 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
538 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
539 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
540 <               rij, pha,dpha,d2pha)
530 >
531 >          ! Calculate rho and drho for atom1
532 >
533 >          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%rho, &
534 >               rij, rha, drha)
535 >          
536 >          ! Calculate Phi(r) for atom1.
537 >          
538 >          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%phi, &
539 >               rij, pha, dpha)
540 >
541         endif
542  
543         if (rij.lt.rcj) then      
544 <          ! Calculate rho,drho and d2rho for atom1
545 <          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
546 <               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
547 <               EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
548 <               EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
549 <               rij, rhb,drhb,d2rhb)
550 <          
551 <          !! Calculate Phi(r) for atom2.
552 <          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
553 <               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
554 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
735 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
736 <               rij, phb,dphb,d2phb)
544 >
545 >          ! Calculate rho and drho for atom2
546 >
547 >          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%rho, &
548 >               rij, rhb, drhb)
549 >
550 >          ! Calculate Phi(r) for atom2.
551 >
552 >          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%phi, &
553 >               rij, phb, dphb)
554 >
555         endif
556  
557         if (rij.lt.rci) then
558            phab = phab + 0.5E0_DP*(rhb/rha)*pha
559            dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
560                 pha*((drhb/rha) - (rhb*drha/rha/rha)))
743 –          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rhb/rha)*d2pha + &
744 –               2.0E0_DP*dpha*((drhb/rha) - (rhb*drha/rha/rha)) + &
745 –               pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
746 –               (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
561         endif
562 <      
562 >
563         if (rij.lt.rcj) then
564            phab = phab + 0.5E0_DP*(rha/rhb)*phb
565            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
566                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
753 –          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rha/rhb)*d2phb + &
754 –               2.0E0_DP*dphb*((drha/rhb) - (rha*drhb/rhb/rhb)) + &
755 –               phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
756 –               (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
567         endif
568 <      
568 >
569         drhoidr = drha
570         drhojdr = drhb
571  
762 –       d2rhoidrdr = d2rha
763 –       d2rhojdrdr = d2rhb
764 –
765 –
572   #ifdef IS_MPI
573         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
574              + dvpdr
# Line 770 | Line 576 | contains
576   #else
577         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
578              + dvpdr
773 –      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
579   #endif
580  
581         fx = dudr * drdx
# Line 780 | Line 585 | contains
585  
586   #ifdef IS_MPI
587         if (do_pot) then
588 <          pot_Row(atom1) = pot_Row(atom1) + phab*0.5
589 <          pot_Col(atom2) = pot_Col(atom2) + phab*0.5
588 >          pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + phab*0.5
589 >          pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + phab*0.5
590         end if
591  
592         f_Row(1,atom1) = f_Row(1,atom1) + fx
593         f_Row(2,atom1) = f_Row(2,atom1) + fy
594         f_Row(3,atom1) = f_Row(3,atom1) + fz
595 <      
595 >
596         f_Col(1,atom2) = f_Col(1,atom2) - fx
597         f_Col(2,atom2) = f_Col(2,atom2) - fy
598         f_Col(3,atom2) = f_Col(3,atom2) - fz
# Line 800 | Line 605 | contains
605         f(1,atom1) = f(1,atom1) + fx
606         f(2,atom1) = f(2,atom1) + fy
607         f(3,atom1) = f(3,atom1) + fz
608 <      
608 >
609         f(1,atom2) = f(1,atom2) - fx
610         f(2,atom2) = f(2,atom2) - fy
611         f(3,atom2) = f(3,atom2) - fz
612   #endif
613 <      
613 >
614         vpair = vpair + phab
615 +
616   #ifdef IS_MPI
617         id1 = AtomRowToGlobal(atom1)
618         id2 = AtomColToGlobal(atom2)
# Line 814 | Line 620 | contains
620         id1 = atom1
621         id2 = atom2
622   #endif
623 <      
623 >
624         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
625 <          
625 >
626            fpair(1) = fpair(1) + fx
627            fpair(2) = fpair(2) + fy
628            fpair(3) = fpair(3) + fz
823 –          
824 –       endif
825 –    
826 –       if (nmflag) then
629  
828 –          drhoidr = drha
829 –          drhojdr = drhb
830 –          d2rhoidrdr = d2rha
831 –          d2rhojdrdr = d2rhb
832 –
833 – #ifdef IS_MPI
834 –          d2 = d2vpdrdr + &
835 –               d2rhoidrdr*dfrhodrho_col(atom2) + &
836 –               d2rhojdrdr*dfrhodrho_row(atom1) + &
837 –               drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
838 –               drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
839 –              
840 – #else
841 –
842 –          d2 = d2vpdrdr + &
843 –               d2rhoidrdr*dfrhodrho(atom2) + &
844 –               d2rhojdrdr*dfrhodrho(atom1) + &
845 –               drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
846 –               drhojdr*drhojdr*d2frhodrhodrho(atom1)
847 – #endif
848 –       end if
849 –      
850 –    endif    
851 –  end subroutine do_eam_pair
852 –
853 –
854 –  subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
855 –
856 –    integer :: atype, nx, j
857 –    real( kind = DP ), dimension(:) :: xa
858 –    real( kind = DP ), dimension(:) :: ya
859 –    real( kind = DP ), dimension(:) :: yppa
860 –    real( kind = DP ) :: x, y
861 –    real( kind = DP ) :: dy, d2y
862 –    real( kind = DP ) :: del, h, a, b, c, d
863 –    integer :: pp_arraySize
864 –
865 –
866 –    ! this spline code assumes that the x points are equally spaced
867 –    ! do not attempt to use this code if they are not.
868 –    
869 –    
870 –    ! find the closest point with a value below our own:
871 –    j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
872 –
873 –    ! check to make sure we're inside the spline range:
874 –    if ((j.gt.nx).or.(j.lt.1)) then
875 –       write(errMSG,*) "EAM_splint: x is outside bounds of spline: ",x,j
876 –       call handleError(routineName,errMSG)
877 –    endif
878 –    ! check to make sure we haven't screwed up the calculation of j:
879 –    if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
880 –       if (j.ne.nx) then
881 –        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
882 –       call handleError(routineName,errMSG)
630         endif
631      endif
632 +  end subroutine do_eam_pair
633  
634 <    del = xa(j+1) - x
887 <    h = xa(j+1) - xa(j)
634 >  subroutine lookupEAMSpline(cs, xval, yval)
635      
636 <    a = del / h
890 <    b = 1.0E0_DP - a
891 <    c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
892 <    d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
893 <    
894 <    y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
895 <  
896 <       dy = (ya(j+1)-ya(j))/h &
897 <            - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
898 <            + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
899 <  
900 <  
901 <       d2y = a*yppa(j) + b*yppa(j+1)
902 <  
636 >    implicit none
637  
638 <  end subroutine eam_splint
639 <
640 <
641 <  subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
642 <
643 <
644 <    ! yp1 and ypn are the first derivatives of y at the two endpoints
645 <    ! if boundary is 'L' the lower derivative is used
912 <    ! if boundary is 'U' the upper derivative is used
913 <    ! if boundary is 'B' then both derivatives are used
914 <    ! if boundary is anything else, then both derivatives are assumed to be 0
638 >    type (cubicSpline), intent(in) :: cs
639 >    real( kind = DP ), intent(in)  :: xval
640 >    real( kind = DP ), intent(out) :: yval
641 >    real( kind = DP ) ::  dx
642 >    integer :: i, j
643 >    !
644 >    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
645 >    !  or is nearest to xval.
646      
647 <    integer :: nx, i, k, max_array_size
647 >    j = MAX(1, MIN(cs%n-1, idint((xval-cs%x(1)) * cs%dx_i) + 1))
648      
649 <    real( kind = DP ), dimension(:)        :: xa
650 <    real( kind = DP ), dimension(:)        :: ya
920 <    real( kind = DP ), dimension(:)        :: yppa
921 <    real( kind = DP ), dimension(size(xa)) :: u
922 <    real( kind = DP ) :: yp1,ypn,un,qn,sig,p
923 <    character(len=*) :: boundary
649 >    dx = xval - cs%x(j)
650 >    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
651      
652 <    ! make sure the sizes match
653 <    if ((nx /= size(xa)) .or. (nx /= size(ya))) then
927 <       call handleWarning("EAM_SPLINE","Array size mismatch")
928 <    end if
652 >    return
653 >  end subroutine lookupEAMSpline
654  
655 <    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
931 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
932 <       yppa(1) = -0.5E0_DP
933 <       u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
934 <            ya(1))/(xa(2)-xa(1))-yp1)
935 <    else
936 <       yppa(1) = 0.0E0_DP
937 <       u(1)  = 0.0E0_DP
938 <    endif
655 >  subroutine lookupEAMSpline1d(cs, xval, yval, dydx)
656      
657 <    do i = 2, nx - 1
941 <       sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
942 <       p = sig * yppa(i-1) + 2.0E0_DP
943 <       yppa(i) = (sig - 1.0E0_DP) / p
944 <       u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
945 <            (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
946 <            (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
947 <    enddo
948 <    
949 <    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
950 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
951 <       qn = 0.5E0_DP
952 <       un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
953 <            (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
954 <    else
955 <       qn = 0.0E0_DP
956 <       un = 0.0E0_DP
957 <    endif
657 >    implicit none
658  
659 <    yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
659 >    type (cubicSpline), intent(in) :: cs
660 >    real( kind = DP ), intent(in)  :: xval
661 >    real( kind = DP ), intent(out) :: yval, dydx
662 >    real( kind = DP ) :: dx
663 >    integer :: i, j
664      
665 <    do k = nx-1, 1, -1
666 <       yppa(k)=yppa(k)*yppa(k+1)+u(k)
963 <    enddo
665 >    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
666 >    !  or is nearest to xval.
667  
965 –  end subroutine eam_spline
668  
669 +    j = MAX(1, MIN(cs%n-1, idint((xval-cs%x(1)) * cs%dx_i) + 1))
670 +    
671 +    dx = xval - cs%x(j)
672 +    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
673  
674 +    dydx = cs%b(j) + dx*(2.0d0 * cs%c(j) + 3.0d0 * dx * cs%d(j))
675 +      
676 +    return
677 +  end subroutine lookupEAMSpline1d
678  
969 –
679   end module eam
971 –
972 –  subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
973 –       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
974 –       eam_ident,status)
975 –       use definitions, ONLY : dp
976 –       use eam, ONLY : module_newEAMtype => newEAMtype
977 –    real (kind = dp )                      :: lattice_constant
978 –    integer                                :: eam_nrho
979 –    real (kind = dp )                      :: eam_drho
980 –    integer                                :: eam_nr
981 –    real (kind = dp )                     :: eam_dr
982 –    real (kind = dp )                      :: rcut
983 –    real (kind = dp ), dimension(eam_nr)    :: eam_Z_r
984 –    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
985 –    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
986 –    integer                                :: eam_ident
987 –    integer                                :: status
988 –
989 –    
990 –    call module_newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
991 –       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
992 –       eam_ident,status)
993 –      
994 – end subroutine newEAMtype

Diff Legend

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