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 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

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

Diff Legend

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