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 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2722 by gezelter, Thu Apr 20 18:24:24 2006 UTC

# Line 44 | Line 44 | module eam
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 65 | Line 68 | module eam
68    character(len=*), parameter :: RoutineName =  "EAM MODULE"
69    !! Logical that determines if eam arrays should be zeroed
70    logical :: cleanme = .true.
68 –  logical :: nmflag  = .false.
71  
70 –
72    type, private :: EAMtype
73       integer           :: eam_atype      
73 –     real( kind = DP ) :: eam_dr          
74 –     integer           :: eam_nr          
75 –     integer           :: eam_nrho          
74       real( kind = DP ) :: eam_lattice        
77 –     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()
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()
77 >     type(cubicSpline) :: rho
78 >     type(cubicSpline) :: Z
79 >     type(cubicSpline) :: F
80 >     type(cubicSpline) :: phi
81    end type EAMtype
82  
93 –
83    !! Arrays for derivatives used in force calculation
84    real( kind = dp), dimension(:), allocatable :: frho
85    real( kind = dp), dimension(:), allocatable :: rho
97 –
86    real( kind = dp), dimension(:), allocatable :: dfrhodrho
99 –  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
87  
101 –
88    !! Arrays for MPI storage
89   #ifdef IS_MPI
90    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
# Line 108 | 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
111 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
112 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
97   #endif
98  
99    type, private :: EAMTypeList
# Line 117 | Line 101 | module eam
101       integer           :: currentAddition = 0
102  
103       type (EAMtype), pointer  :: EAMParams(:) => null()
104 +     integer, pointer         :: atidToEAMType(:) => null()
105    end type EAMTypeList
106  
122 –
107    type (eamTypeList), save :: EAMList
108  
109    !! standard eam stuff  
# Line 133 | Line 117 | module eam
117    public :: calc_eam_preforce_Frho
118    public :: clean_EAM
119    public :: destroyEAMTypes
120 +  public :: getEAMCut
121  
122   contains
123  
139 –
124    subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
125 <       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
142 <       eam_ident,status)
125 >       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho, c_ident, status)
126      real (kind = dp )                      :: lattice_constant
127      integer                                :: eam_nrho
128      real (kind = dp )                      :: eam_drho
129      integer                                :: eam_nr
130      real (kind = dp )                      :: eam_dr
131      real (kind = dp )                      :: rcut
132 <    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r
133 <    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
134 <    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
135 <    integer                                :: eam_ident
132 >    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r, rvals
133 >    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r, eam_phi_r
134 >    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho, rhovals
135 >    integer                                :: c_ident
136      integer                                :: status
137  
138 <    integer                                :: nAtypes
138 >    integer                                :: nAtypes,nEAMTypes,myATID
139      integer                                :: maxVals
140      integer                                :: alloc_stat
141 <    integer                                :: current
141 >    integer                                :: current, j
142      integer,pointer                        :: Matchlist(:) => null()
143  
144      status = 0
145  
163 –
146      !! Assume that atypes has already been set and get the total number of types in atypes
147      !! Also assume that every member of atypes is a EAM model.
148  
167 –
149      ! check to see if this is the first time into
150      if (.not.associated(EAMList%EAMParams)) then
151 <       call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
152 <       EAMList%n_eam_types = nAtypes
153 <       allocate(EAMList%EAMParams(nAtypes))
151 >       call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
152 >       EAMList%n_eam_types = nEAMtypes
153 >       allocate(EAMList%EAMParams(nEAMTypes))
154 >       nAtypes = getSize(atypes)
155 >       allocate(EAMList%atidToEAMType(nAtypes))
156      end if
157  
158      EAMList%currentAddition = EAMList%currentAddition + 1
159      current = EAMList%currentAddition
160  
161 +    myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
162 +    EAMList%atidToEAMType(myATID) = current
163  
164 <    call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
180 <    if (alloc_stat /= 0) then
181 <       status = -1
182 <       return
183 <    end if
184 <
185 <    ! this is a possible bug, we assume a correspondence between the vector atypes and
186 <    ! EAMAtypes
187 <
188 <    EAMList%EAMParams(current)%eam_atype    = eam_ident
164 >    EAMList%EAMParams(current)%eam_atype    = c_ident
165      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
190 –    EAMList%EAMParams(current)%eam_nrho     = eam_nrho
191 –    EAMList%EAMParams(current)%eam_drho     = eam_drho
192 –    EAMList%EAMParams(current)%eam_nr       = eam_nr
193 –    EAMList%EAMParams(current)%eam_dr       = eam_dr
166      EAMList%EAMParams(current)%eam_rcut     = rcut
167 <    EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
168 <    EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
169 <    EAMList%EAMParams(current)%eam_F_rho    = eam_F_rho
167 >
168 >    ! Build array of r values
169 >    do j = 1, eam_nr
170 >       rvals(j) = real(j-1,kind=dp) * eam_dr
171 >    end do
172 >    
173 >    ! Build array of rho values
174 >    do j = 1, eam_nrho
175 >       rhovals(j) = real(j-1,kind=dp) * eam_drho
176 >    end do
177  
178 +    ! convert from eV to kcal / mol:
179 +    do j = 1, eam_nrho
180 +       eam_F_rho(j) = eam_F_rho(j) * 23.06054E0_DP
181 +    end do
182 +    
183 +    ! precompute the pair potential and get it into kcal / mol:
184 +    eam_phi_r(1) = 0.0E0_DP
185 +    do j = 2, eam_nr
186 +       eam_phi_r(j) = 331.999296E0_DP * (eam_Z_r(j)**2) / rvals(j)
187 +    enddo
188 +
189 +    call newSpline(EAMList%EAMParams(current)%rho, rvals, eam_rho_r, .true.)
190 +    call newSpline(EAMList%EAMParams(current)%Z,   rvals, eam_Z_r, .true.)
191 +    call newSpline(EAMList%EAMParams(current)%F, rhovals, eam_F_rho, .true.)
192 +    call newSpline(EAMList%EAMParams(current)%phi, rvals, eam_phi_r, .true.)
193    end subroutine newEAMtype
194  
195  
# Line 213 | Line 207 | contains
207  
208      eamList%n_eam_types = 0
209      eamList%currentAddition = 0
216 –
210    end subroutine destroyEAMtypes
211  
212 +  function getEAMCut(atomID) result(cutValue)
213 +    integer, intent(in) :: atomID
214 +    integer :: eamID
215 +    real(kind=dp) :: cutValue
216 +    
217 +    eamID = EAMList%atidToEAMType(atomID)
218 +    cutValue = EAMList%EAMParams(eamID)%eam_rcut
219 +  end function getEAMCut
220 +
221    subroutine init_EAM_FF(status)
222      integer :: status
223      integer :: i,j
224      real(kind=dp) :: current_rcut_max
225 + #ifdef IS_MPI
226 +    integer :: nAtomsInRow
227 +    integer :: nAtomsInCol
228 + #endif
229      integer :: alloc_stat
230      integer :: number_r, number_rho
231  
226 –
232      status = 0
233      if (EAMList%currentAddition == 0) then
234         call handleError("init_EAM_FF","No members in EAMList")
235         status = -1
236         return
237      end if
238 <
234 <
235 <    do i = 1, EAMList%currentAddition
236 <
237 <       ! Build array of r values
238 <
239 <       do j = 1,EAMList%EAMParams(i)%eam_nr
240 <          EAMList%EAMParams(i)%eam_rvals(j) = &
241 <               real(j-1,kind=dp)* &
242 <               EAMList%EAMParams(i)%eam_dr
243 <       end do
244 <       ! Build array of rho values
245 <       do j = 1,EAMList%EAMParams(i)%eam_nrho
246 <          EAMList%EAMParams(i)%eam_rhovals(j) = &
247 <               real(j-1,kind=dp)* &
248 <               EAMList%EAMParams(i)%eam_drho
249 <       end do
250 <       ! convert from eV to kcal / mol:
251 <       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
252 <
253 <       ! precompute the pair potential and get it into kcal / mol:
254 <       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
255 <       do j = 2, EAMList%EAMParams(i)%eam_nr
256 <          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
257 <          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
258 <       enddo
259 <    end do
260 <
261 <
262 <    do i = 1,  EAMList%currentAddition
263 <       number_r   = EAMList%EAMParams(i)%eam_nr
264 <       number_rho = EAMList%EAMParams(i)%eam_nrho
265 <
266 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
267 <            EAMList%EAMParams(i)%eam_rho_r, &
268 <            EAMList%EAMParams(i)%eam_rho_r_pp, &
269 <            0.0E0_DP, 0.0E0_DP, 'N')
270 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
271 <            EAMList%EAMParams(i)%eam_Z_r, &
272 <            EAMList%EAMParams(i)%eam_Z_r_pp, &
273 <            0.0E0_DP, 0.0E0_DP, 'N')
274 <       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
275 <            EAMList%EAMParams(i)%eam_F_rho, &
276 <            EAMList%EAMParams(i)%eam_F_rho_pp, &
277 <            0.0E0_DP, 0.0E0_DP, 'N')
278 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
279 <            EAMList%EAMParams(i)%eam_phi_r, &
280 <            EAMList%EAMParams(i)%eam_phi_r_pp, &
281 <            0.0E0_DP, 0.0E0_DP, 'N')
282 <    enddo
283 <
284 <    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
285 <    !! find the smallest rcut for any eam atype
286 <    !       do i = 2, EAMList%currentAddition
287 <    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
288 <    !       end do
289 <
290 <    !       EAM_rcut = current_rcut_max
291 <    !       EAM_rcut_orig = current_rcut_max
292 <    !       do i = 1, EAMList%currentAddition
293 <    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
294 <    !       end do
238 >    
239      !! Allocate arrays for force calculation
240 <
297 <    call allocateEAM(alloc_stat)
298 <    if (alloc_stat /= 0 ) then
299 <       write(*,*) "allocateEAM failed"
300 <       status = -1
301 <       return
302 <    endif
303 <
304 <  end subroutine init_EAM_FF
305 <
306 <  !! routine checks to see if array is allocated, deallocates array if allocated
307 <  !! and then creates the array to the required size
308 <  subroutine allocateEAM(status)
309 <    integer, intent(out) :: status
310 <
240 >    
241   #ifdef IS_MPI
312 –    integer :: nAtomsInRow
313 –    integer :: nAtomsInCol
314 – #endif
315 –    integer :: alloc_stat
316 –
317 –
318 –    status = 0
319 – #ifdef IS_MPI
242      nAtomsInRow = getNatomsInRow(plan_atom_row)
243      nAtomsInCol = getNatomsInCol(plan_atom_col)
244   #endif
# Line 327 | Line 249 | contains
249         status = -1
250         return
251      end if
252 +
253      if (allocated(rho)) deallocate(rho)
254      allocate(rho(nlocal),stat=alloc_stat)
255      if (alloc_stat /= 0) then
# Line 341 | Line 264 | contains
264         return
265      end if
266  
344 –    if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
345 –    allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
346 –    if (alloc_stat /= 0) then
347 –       status = -1
348 –       return
349 –    end if
350 –
267   #ifdef IS_MPI
268  
269      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 357 | Line 273 | contains
273         return
274      end if
275  
360 –
276      if (allocated(frho_row)) deallocate(frho_row)
277      allocate(frho_row(nAtomsInRow),stat=alloc_stat)
278      if (alloc_stat /= 0) then
# Line 376 | Line 291 | contains
291         status = -1
292         return
293      end if
379 –    if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
380 –    allocate(d2frhodrhodrho_row(nAtomsInRow),stat=alloc_stat)
381 –    if (alloc_stat /= 0) then
382 –       status = -1
383 –       return
384 –    end if
294  
386 –
295      ! Now do column arrays
296  
297      if (allocated(frho_col)) deallocate(frho_col)
# Line 404 | Line 312 | contains
312         status = -1
313         return
314      end if
407 –    if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
408 –    allocate(d2frhodrhodrho_col(nAtomsInCol),stat=alloc_stat)
409 –    if (alloc_stat /= 0) then
410 –       status = -1
411 –       return
412 –    end if
315  
316   #endif
317  
318 <  end subroutine allocateEAM
318 >  end subroutine init_EAM_FF
319  
320 <  !! C sets rcut to be the largest cutoff of any atype
419 <  !! present in this simulation. Doesn't include all atypes
420 <  !! sim knows about, just those in the simulation.
421 <  subroutine setCutoffEAM(rcut, status)
320 >  subroutine setCutoffEAM(rcut)
321      real(kind=dp) :: rcut
423 –    integer :: status
424 –    status = 0
425 –
322      EAM_rcut = rcut
427 –
323    end subroutine setCutoffEAM
324  
430 –
431 –
325    subroutine clean_EAM()
326  
327      ! clean non-IS_MPI first
# Line 447 | Line 340 | contains
340   #endif
341    end subroutine clean_EAM
342  
450 –
451 –
452 –  subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
453 –    integer, intent(in)          :: eam_n_rho
454 –    integer, intent(in)          :: eam_n_r
455 –    type (EAMType)               :: thisEAMType
456 –    integer, optional   :: stat
457 –    integer             :: alloc_stat
458 –
459 –
460 –
461 –    if (present(stat)) stat = 0
462 –
463 –    allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
464 –    if (alloc_stat /= 0 ) then
465 –       if (present(stat)) stat = -1
466 –       return
467 –    end if
468 –    allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat)  
469 –    if (alloc_stat /= 0 ) then
470 –       if (present(stat)) stat = -1
471 –       return
472 –    end if
473 –    allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)  
474 –    if (alloc_stat /= 0 ) then
475 –       if (present(stat)) stat = -1
476 –       return
477 –    end if
478 –    allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)        
479 –    if (alloc_stat /= 0 ) then
480 –       if (present(stat)) stat = -1
481 –       return
482 –    end if
483 –    allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat)      
484 –    if (alloc_stat /= 0 ) then
485 –       if (present(stat)) stat = -1
486 –       return
487 –    end if
488 –    allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat)      
489 –    if (alloc_stat /= 0 ) then
490 –       if (present(stat)) stat = -1
491 –       return
492 –    end if
493 –    allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat)  
494 –    if (alloc_stat /= 0 ) then
495 –       if (present(stat)) stat = -1
496 –       return
497 –    end if
498 –    allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)  
499 –    if (alloc_stat /= 0 ) then
500 –       if (present(stat)) stat = -1
501 –       return
502 –    end if
503 –    allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)  
504 –    if (alloc_stat /= 0 ) then
505 –       if (present(stat)) stat = -1
506 –       return
507 –    end if
508 –    allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
509 –    if (alloc_stat /= 0 ) then
510 –       if (present(stat)) stat = -1
511 –       return
512 –    end if
513 –
514 –
515 –  end subroutine allocate_EAMType
516 –
517 –
343    subroutine deallocate_EAMType(thisEAMType)
344      type (EAMtype), pointer :: thisEAMType
345  
346 <    ! free Arrays in reverse order of allocation...
347 <    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
348 <    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
349 <    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
525 <    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
526 <    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
527 <    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
528 <    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
529 <    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
530 <    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
531 <    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
346 >    call deleteSpline(thisEAMType%F)
347 >    call deleteSpline(thisEAMType%rho)
348 >    call deleteSpline(thisEAMType%phi)
349 >    call deleteSpline(thisEAMType%Z)
350  
351    end subroutine deallocate_EAMType
352  
353    !! Calculates rho_r
354    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
355 <    integer :: atom1,atom2
355 >    integer :: atom1, atom2
356      real(kind = dp), dimension(3) :: d
357      real(kind = dp), intent(inout)               :: r
358      real(kind = dp), intent(inout)               :: rijsq
# Line 542 | Line 360 | contains
360      real(kind = dp) :: rho_i_at_j
361      ! value of electron density rho do to atom j at atom i
362      real(kind = dp) :: rho_j_at_i
545 –
546 –    ! we don't use the derivatives, dummy variables
547 –    real( kind = dp) :: drho,d2rho
363      integer :: eam_err
364 +    
365 +    integer :: atid1, atid2 ! Global atid    
366 +    integer :: myid_atom1 ! EAM atid
367 +    integer :: myid_atom2
368  
550 –    integer :: myid_atom1
551 –    integer :: myid_atom2
552 –
369      ! check to see if we need to be cleaned at the start of a force loop
370  
555 –
556 –
557 –
371   #ifdef IS_MPI
372 <    myid_atom1 = atid_Row(atom1)
373 <    myid_atom2 = atid_Col(atom2)
372 >    Atid1 = Atid_row(Atom1)
373 >    Atid2 = Atid_col(Atom2)
374   #else
375 <    myid_atom1 = atid(atom1)
376 <    myid_atom2 = atid(atom2)
375 >    Atid1 = Atid(Atom1)
376 >    Atid2 = Atid(Atom2)
377   #endif
378  
379 +    Myid_atom1 = Eamlist%atidtoeamtype(Atid1)
380 +    Myid_atom2 = Eamlist%atidtoeamtype(Atid2)
381 +
382      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
383  
384 +       call lookupUniformSpline(EAMList%EAMParams(myid_atom1)%rho, r, &
385 +            rho_i_at_j)
386  
569 –
570 –       call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
571 –            EAMList%EAMParams(myid_atom1)%eam_rvals, &
572 –            EAMList%EAMParams(myid_atom1)%eam_rho_r, &
573 –            EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
574 –            r, rho_i_at_j,drho,d2rho)
575 –
576 –
577 –
387   #ifdef  IS_MPI
388         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
389   #else
390         rho(atom2) = rho(atom2) + rho_i_at_j
391   #endif
583 –       !       write(*,*) atom1,atom2,r,rho_i_at_j
392      endif
393  
394      if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
587 –       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
588 –            EAMList%EAMParams(myid_atom2)%eam_rvals, &
589 –            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
590 –            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
591 –            r, rho_j_at_i,drho,d2rho)
395  
396 +       call lookupUniformSpline(EAMList%EAMParams(myid_atom2)%rho, r, &
397 +            rho_j_at_i)
398  
594 –
595 –
399   #ifdef  IS_MPI
400         rho_row(atom1) = rho_row(atom1) + rho_j_at_i
401   #else
402         rho(atom1) = rho(atom1) + rho_j_at_i
403   #endif
404      endif
602 –
603 –
604 –
605 –
606 –
607 –
405    end subroutine calc_eam_prepair_rho
406  
407  
611 –
612 –
408    !! Calculate the functional F(rho) for all local atoms
409 <  subroutine calc_eam_preforce_Frho(nlocal,pot)
409 >  subroutine calc_eam_preforce_Frho(nlocal, pot)
410      integer :: nlocal
411      real(kind=dp) :: pot
412 <    integer :: i,j
412 >    integer :: i, j
413      integer :: atom
414 <    real(kind=dp) :: U,U1,U2
414 >    real(kind=dp) :: U,U1
415      integer :: atype1
416 <    integer :: me
622 <    integer :: n_rho_points
416 >    integer :: me, atid1
417  
624 –
418      cleanme = .true.
419 <    !! Scatter the electron density from  pre-pair calculation back to local atoms
419 >    !! Scatter the electron density from  pre-pair calculation back to
420 >    !! local atoms
421   #ifdef IS_MPI
422      call scatter(rho_row,rho,plan_atom_row,eam_err)
423      if (eam_err /= 0 ) then
# Line 639 | Line 433 | contains
433      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
434   #endif
435  
642 –
643 –
436      !! Calculate F(rho) and derivative
437      do atom = 1, nlocal
438 <       me = atid(atom)
439 <       n_rho_points = EAMList%EAMParams(me)%eam_nrho
648 <       !  Check to see that the density is not greater than the larges rho we have calculated
649 <       if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
650 <          call eam_splint(n_rho_points, &
651 <               EAMList%EAMParams(me)%eam_rhovals, &
652 <               EAMList%EAMParams(me)%eam_f_rho, &
653 <               EAMList%EAMParams(me)%eam_f_rho_pp, &
654 <               rho(atom), & ! Actual Rho
655 <               u, u1, u2)
656 <       else
657 <          ! Calculate F(rho with the largest available rho value
658 <          call eam_splint(n_rho_points, &
659 <               EAMList%EAMParams(me)%eam_rhovals, &
660 <               EAMList%EAMParams(me)%eam_f_rho, &
661 <               EAMList%EAMParams(me)%eam_f_rho_pp, &
662 <               EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
663 <               u,u1,u2)
664 <       end if
438 >       atid1 = atid(atom)
439 >       me = eamList%atidToEAMtype(atid1)
440  
441 <
441 >       call lookupUniformSpline1d(EAMList%EAMParams(me)%F, rho(atom), &
442 >            u, u1)
443 >      
444         frho(atom) = u
445         dfrhodrho(atom) = u1
669 –       d2frhodrhodrho(atom) = u2
446         pot = pot + u
447  
448      enddo
449  
674 –
675 –
450   #ifdef IS_MPI
451      !! communicate f(rho) and derivatives back into row and column arrays
452      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 691 | Line 465 | contains
465      if (eam_err /=  0) then
466         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
467      endif
694 –
695 –
696 –
697 –
698 –
699 –    if (nmflag) then
700 –       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row)
701 –       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col)
702 –    endif
468   #endif
469  
470  
471    end subroutine calc_eam_preforce_Frho
472 <
708 <
709 <
710 <
472 >  
473    !! Does EAM pairwise Force calculation.  
474    subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
475         pot, f, do_pot)
# Line 721 | Line 483 | contains
483  
484      logical, intent(in) :: do_pot
485  
486 <    real( kind = dp ) :: drdx,drdy,drdz
487 <    real( kind = dp ) :: d2
488 <    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
489 <    real( kind = dp ) :: rha,drha,d2rha, dpha
728 <    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
486 >    real( kind = dp ) :: drdx, drdy, drdz
487 >    real( kind = dp ) :: phab, pha, dvpdr
488 >    real( kind = dp ) :: rha, drha, dpha
489 >    real( kind = dp ) :: rhb, drhb, dphb
490      real( kind = dp ) :: dudr
491 <    real( kind = dp ) :: rci,rcj
492 <    real( kind = dp ) :: drhoidr,drhojdr
493 <    real( kind = dp ) :: d2rhoidrdr
494 <    real( kind = dp ) :: d2rhojdrdr
734 <    real( kind = dp ) :: Fx,Fy,Fz
735 <    real( kind = dp ) :: r,d2pha,phb,d2phb
491 >    real( kind = dp ) :: rci, rcj
492 >    real( kind = dp ) :: drhoidr, drhojdr
493 >    real( kind = dp ) :: Fx, Fy, Fz
494 >    real( kind = dp ) :: r, phb
495  
496 <    integer :: id1,id2
496 >    integer :: id1, id2
497      integer  :: mytype_atom1
498      integer  :: mytype_atom2
499 +    integer  :: atid1, atid2
500  
741 –    !Local Variables
742 –
743 –    ! write(*,*) "Frho: ", Frho(atom1)
744 –
501      phab = 0.0E0_DP
502      dvpdr = 0.0E0_DP
747 –    d2vpdrdr = 0.0E0_DP
503  
504      if (rij .lt. EAM_rcut) then
505  
506   #ifdef IS_MPI
507 <       mytype_atom1 = atid_row(atom1)
508 <       mytype_atom2 = atid_col(atom2)
507 >       atid1 = atid_row(atom1)
508 >       atid2 = atid_col(atom2)
509   #else
510 <       mytype_atom1 = atid(atom1)
511 <       mytype_atom2 = atid(atom2)
510 >       atid1 = atid(atom1)
511 >       atid2 = atid(atom2)
512   #endif
513 +
514 +       mytype_atom1 = EAMList%atidToEAMType(atid1)
515 +       mytype_atom2 = EAMList%atidTOEAMType(atid2)
516 +
517 +
518         ! get cutoff for atom 1
519         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
520         ! get type specific cutoff for atom 2
# Line 765 | Line 525 | contains
525         drdz = d(3)/rij
526  
527         if (rij.lt.rci) then
528 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
529 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
530 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
531 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
532 <               rij, rha,drha,d2rha)
533 <          !! Calculate Phi(r) for atom1.
534 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
535 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
536 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
537 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
538 <               rij, pha,dpha,d2pha)
528 >
529 >          ! Calculate rho and drho for atom1
530 >
531 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom1)%rho, &
532 >               rij, rha, drha)
533 >          
534 >          ! Calculate Phi(r) for atom1.
535 >          
536 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom1)%phi, &
537 >               rij, pha, dpha)
538 >
539         endif
540  
541         if (rij.lt.rcj) then      
782 –          ! Calculate rho,drho and d2rho for atom1
783 –          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
784 –               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
785 –               EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
786 –               EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
787 –               rij, rhb,drhb,d2rhb)
542  
543 <          !! Calculate Phi(r) for atom2.
544 <          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
545 <               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
546 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
547 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
548 <               rij, phb,dphb,d2phb)
543 >          ! Calculate rho and drho for atom2
544 >
545 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom2)%rho, &
546 >               rij, rhb, drhb)
547 >
548 >          ! Calculate Phi(r) for atom2.
549 >
550 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom2)%phi, &
551 >               rij, phb, dphb)
552 >
553         endif
554  
555         if (rij.lt.rci) then
556            phab = phab + 0.5E0_DP*(rhb/rha)*pha
557            dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
558                 pha*((drhb/rha) - (rhb*drha/rha/rha)))
801 –          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rhb/rha)*d2pha + &
802 –               2.0E0_DP*dpha*((drhb/rha) - (rhb*drha/rha/rha)) + &
803 –               pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
804 –               (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
559         endif
560  
561         if (rij.lt.rcj) then
562            phab = phab + 0.5E0_DP*(rha/rhb)*phb
563            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
564                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
811 –          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rha/rhb)*d2phb + &
812 –               2.0E0_DP*dphb*((drha/rhb) - (rha*drhb/rhb/rhb)) + &
813 –               phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
814 –               (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
565         endif
566  
567         drhoidr = drha
568         drhojdr = drhb
569  
820 –       d2rhoidrdr = d2rha
821 –       d2rhojdrdr = d2rhb
822 –
823 –
570   #ifdef IS_MPI
571         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
572              + dvpdr
# Line 828 | Line 574 | contains
574   #else
575         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
576              + dvpdr
831 –       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
577   #endif
578  
579         fx = dudr * drdx
# Line 838 | Line 583 | contains
583  
584   #ifdef IS_MPI
585         if (do_pot) then
586 <          pot_Row(atom1) = pot_Row(atom1) + phab*0.5
587 <          pot_Col(atom2) = pot_Col(atom2) + phab*0.5
586 >          pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + phab*0.5
587 >          pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + phab*0.5
588         end if
589  
590         f_Row(1,atom1) = f_Row(1,atom1) + fx
# Line 865 | Line 610 | contains
610   #endif
611  
612         vpair = vpair + phab
613 +
614   #ifdef IS_MPI
615         id1 = AtomRowToGlobal(atom1)
616         id2 = AtomColToGlobal(atom2)
# Line 880 | Line 626 | contains
626            fpair(3) = fpair(3) + fz
627  
628         endif
883 –
884 –       if (nmflag) then
885 –
886 –          drhoidr = drha
887 –          drhojdr = drhb
888 –          d2rhoidrdr = d2rha
889 –          d2rhojdrdr = d2rhb
890 –
891 – #ifdef IS_MPI
892 –          d2 = d2vpdrdr + &
893 –               d2rhoidrdr*dfrhodrho_col(atom2) + &
894 –               d2rhojdrdr*dfrhodrho_row(atom1) + &
895 –               drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
896 –               drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
897 –
898 – #else
899 –
900 –          d2 = d2vpdrdr + &
901 –               d2rhoidrdr*dfrhodrho(atom2) + &
902 –               d2rhojdrdr*dfrhodrho(atom1) + &
903 –               drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
904 –               drhojdr*drhojdr*d2frhodrhodrho(atom1)
905 – #endif
906 –       end if
907 –
629      endif
630    end subroutine do_eam_pair
631  
911 –
912 –  subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
913 –
914 –    integer :: atype, nx, j
915 –    real( kind = DP ), dimension(:) :: xa
916 –    real( kind = DP ), dimension(:) :: ya
917 –    real( kind = DP ), dimension(:) :: yppa
918 –    real( kind = DP ) :: x, y
919 –    real( kind = DP ) :: dy, d2y
920 –    real( kind = DP ) :: del, h, a, b, c, d
921 –    integer :: pp_arraySize
922 –
923 –
924 –    ! this spline code assumes that the x points are equally spaced
925 –    ! do not attempt to use this code if they are not.
926 –
927 –
928 –    ! find the closest point with a value below our own:
929 –    j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
930 –
931 –    ! check to make sure we're inside the spline range:
932 –    if ((j.gt.nx).or.(j.lt.1)) then
933 –       write(errMSG,*) "EAM_splint: x is outside bounds of spline: ",x,j
934 –       call handleError(routineName,errMSG)
935 –    endif
936 –    ! check to make sure we haven't screwed up the calculation of j:
937 –    if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
938 –       if (j.ne.nx) then
939 –          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
940 –          call handleError(routineName,errMSG)
941 –       endif
942 –    endif
943 –
944 –    del = xa(j+1) - x
945 –    h = xa(j+1) - xa(j)
946 –
947 –    a = del / h
948 –    b = 1.0E0_DP - a
949 –    c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
950 –    d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
951 –
952 –    y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
953 –
954 –    dy = (ya(j+1)-ya(j))/h &
955 –         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
956 –         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
957 –
958 –
959 –    d2y = a*yppa(j) + b*yppa(j+1)
960 –
961 –
962 –  end subroutine eam_splint
963 –
964 –
965 –  subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
966 –
967 –
968 –    ! yp1 and ypn are the first derivatives of y at the two endpoints
969 –    ! if boundary is 'L' the lower derivative is used
970 –    ! if boundary is 'U' the upper derivative is used
971 –    ! if boundary is 'B' then both derivatives are used
972 –    ! if boundary is anything else, then both derivatives are assumed to be 0
973 –
974 –    integer :: nx, i, k, max_array_size
975 –
976 –    real( kind = DP ), dimension(:)        :: xa
977 –    real( kind = DP ), dimension(:)        :: ya
978 –    real( kind = DP ), dimension(:)        :: yppa
979 –    real( kind = DP ), dimension(size(xa)) :: u
980 –    real( kind = DP ) :: yp1,ypn,un,qn,sig,p
981 –    character(len=*) :: boundary
982 –
983 –    ! make sure the sizes match
984 –    if ((nx /= size(xa)) .or. (nx /= size(ya))) then
985 –       call handleWarning("EAM_SPLINE","Array size mismatch")
986 –    end if
987 –
988 –    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
989 –         (boundary.eq.'b').or.(boundary.eq.'B')) then
990 –       yppa(1) = -0.5E0_DP
991 –       u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
992 –            ya(1))/(xa(2)-xa(1))-yp1)
993 –    else
994 –       yppa(1) = 0.0E0_DP
995 –       u(1)  = 0.0E0_DP
996 –    endif
997 –
998 –    do i = 2, nx - 1
999 –       sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1000 –       p = sig * yppa(i-1) + 2.0E0_DP
1001 –       yppa(i) = (sig - 1.0E0_DP) / p
1002 –       u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
1003 –            (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1004 –            (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1005 –    enddo
1006 –
1007 –    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1008 –         (boundary.eq.'b').or.(boundary.eq.'B')) then
1009 –       qn = 0.5E0_DP
1010 –       un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
1011 –            (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
1012 –    else
1013 –       qn = 0.0E0_DP
1014 –       un = 0.0E0_DP
1015 –    endif
1016 –
1017 –    yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1018 –
1019 –    do k = nx-1, 1, -1
1020 –       yppa(k)=yppa(k)*yppa(k+1)+u(k)
1021 –    enddo
1022 –
1023 –  end subroutine eam_spline
1024 –
632   end module eam

Diff Legend

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