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 2716 by gezelter, Wed Oct 12 21:00:50 2005 UTC vs.
Revision 2717 by gezelter, Mon Apr 17 21:49:12 2006 UTC

# Line 45 | Line 45 | module eam
45    use status
46    use atype_module
47    use vector_class
48 +  use interpolation
49   #ifdef IS_MPI
50    use mpiSimulation
51   #endif
# Line 67 | 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.
70 –  logical :: nmflag  = .false.
71  
72 –
72    type, private :: EAMtype
73       integer           :: eam_atype      
75 –     real( kind = DP ) :: eam_dr          
76 –     integer           :: eam_nr          
77 –     integer           :: eam_nrho          
74       real( kind = DP ) :: eam_lattice        
79 –     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()
86 <     real( kind = DP ), pointer, dimension(:) :: eam_Z_r          => null()
87 <     real( kind = DP ), pointer, dimension(:) :: eam_rho_r        => null()
88 <     real( kind = DP ), pointer, dimension(:) :: eam_phi_r        => null()
89 <     real( kind = DP ), pointer, dimension(:) :: eam_F_rho_pp     => null()
90 <     real( kind = DP ), pointer, dimension(:) :: eam_Z_r_pp       => null()
91 <     real( kind = DP ), pointer, dimension(:) :: eam_rho_r_pp     => null()
92 <     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  
95 –
83    !! Arrays for derivatives used in force calculation
84    real( kind = dp), dimension(:), allocatable :: frho
85    real( kind = dp), dimension(:), allocatable :: rho
99 –
86    real( kind = dp), dimension(:), allocatable :: dfrhodrho
101 –  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
87  
103 –
88    !! Arrays for MPI storage
89   #ifdef IS_MPI
90    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
# Line 110 | 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
113 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
114 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
97   #endif
98  
99    type, private :: EAMTypeList
# Line 122 | Line 104 | module eam
104       integer, pointer         :: atidToEAMType(:) => null()
105    end type EAMTypeList
106  
125 –
107    type (eamTypeList), save :: EAMList
108  
109    !! standard eam stuff  
# Line 140 | Line 121 | contains
121  
122   contains
123  
143 –
124    subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
125 <       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
146 <       c_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
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,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  
167 –
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  
171 –
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., nEAMtypes, MatchList)
# Line 184 | Line 161 | contains
161      myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
162      EAMList%atidToEAMType(myATID) = current
163  
187 –    call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
188 –    if (alloc_stat /= 0) then
189 –       status = -1
190 –       return
191 –    end if
192 –
193 –  
164      EAMList%EAMParams(current)%eam_atype    = c_ident
165      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
196 –    EAMList%EAMParams(current)%eam_nrho     = eam_nrho
197 –    EAMList%EAMParams(current)%eam_drho     = eam_drho
198 –    EAMList%EAMParams(current)%eam_nr       = eam_nr
199 –    EAMList%EAMParams(current)%eam_dr       = eam_dr
166      EAMList%EAMParams(current)%eam_rcut     = rcut
201 –    EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
202 –    EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
203 –    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, rhovals, &
190 +         0.0d0, 0.0d0, .true.)
191 +    call newSpline(EAMList%EAMParams(current)%Z,   rvals, eam_Z_r, &
192 +         0.0d0, 0.0d0, .true.)
193 +    call newSpline(EAMList%EAMParams(current)%F, rhovals, eam_F_rho, &
194 +         0.0d0, 0.0d0, .true.)
195 +    call newSpline(EAMList%EAMParams(current)%phi, rvals, eam_phi_r, &
196 +         0.0d0, 0.0d0, .true.)
197    end subroutine newEAMtype
198  
199  
# Line 219 | Line 211 | contains
211  
212      eamList%n_eam_types = 0
213      eamList%currentAddition = 0
222 –
214    end subroutine destroyEAMtypes
215  
216    function getEAMCut(atomID) result(cutValue)
# Line 235 | Line 226 | contains
226      integer :: status
227      integer :: i,j
228      real(kind=dp) :: current_rcut_max
229 + #ifdef IS_MPI
230 +    integer :: nAtomsInRow
231 +    integer :: nAtomsInCol
232 + #endif
233      integer :: alloc_stat
234      integer :: number_r, number_rho
235  
241 –
236      status = 0
237      if (EAMList%currentAddition == 0) then
238         call handleError("init_EAM_FF","No members in EAMList")
239         status = -1
240         return
241      end if
242 <
249 <
250 <    do i = 1, EAMList%currentAddition
251 <
252 <       ! Build array of r values
253 <
254 <       do j = 1,EAMList%EAMParams(i)%eam_nr
255 <          EAMList%EAMParams(i)%eam_rvals(j) = &
256 <               real(j-1,kind=dp)* &
257 <               EAMList%EAMParams(i)%eam_dr
258 <       end do
259 <       ! Build array of rho values
260 <       do j = 1,EAMList%EAMParams(i)%eam_nrho
261 <          EAMList%EAMParams(i)%eam_rhovals(j) = &
262 <               real(j-1,kind=dp)* &
263 <               EAMList%EAMParams(i)%eam_drho
264 <       end do
265 <       ! convert from eV to kcal / mol:
266 <       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
267 <
268 <       ! precompute the pair potential and get it into kcal / mol:
269 <       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
270 <       do j = 2, EAMList%EAMParams(i)%eam_nr
271 <          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
272 <          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
273 <       enddo
274 <    end do
275 <
276 <
277 <    do i = 1,  EAMList%currentAddition
278 <       number_r   = EAMList%EAMParams(i)%eam_nr
279 <       number_rho = EAMList%EAMParams(i)%eam_nrho
280 <
281 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
282 <            EAMList%EAMParams(i)%eam_rho_r, &
283 <            EAMList%EAMParams(i)%eam_rho_r_pp, &
284 <            0.0E0_DP, 0.0E0_DP, 'N')
285 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
286 <            EAMList%EAMParams(i)%eam_Z_r, &
287 <            EAMList%EAMParams(i)%eam_Z_r_pp, &
288 <            0.0E0_DP, 0.0E0_DP, 'N')
289 <       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
290 <            EAMList%EAMParams(i)%eam_F_rho, &
291 <            EAMList%EAMParams(i)%eam_F_rho_pp, &
292 <            0.0E0_DP, 0.0E0_DP, 'N')
293 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
294 <            EAMList%EAMParams(i)%eam_phi_r, &
295 <            EAMList%EAMParams(i)%eam_phi_r_pp, &
296 <            0.0E0_DP, 0.0E0_DP, 'N')
297 <    enddo
298 <
299 <    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
300 <    !! find the smallest rcut for any eam atype
301 <    !       do i = 2, EAMList%currentAddition
302 <    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
303 <    !       end do
304 <
305 <    !       EAM_rcut = current_rcut_max
306 <    !       EAM_rcut_orig = current_rcut_max
307 <    !       do i = 1, EAMList%currentAddition
308 <    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
309 <    !       end do
242 >    
243      !! Allocate arrays for force calculation
244 <
312 <    call allocateEAM(alloc_stat)
313 <    if (alloc_stat /= 0 ) then
314 <       write(*,*) "allocateEAM failed"
315 <       status = -1
316 <       return
317 <    endif
318 <
319 <  end subroutine init_EAM_FF
320 <
321 <  !! routine checks to see if array is allocated, deallocates array if allocated
322 <  !! and then creates the array to the required size
323 <  subroutine allocateEAM(status)
324 <    integer, intent(out) :: status
325 <
244 >    
245   #ifdef IS_MPI
327 –    integer :: nAtomsInRow
328 –    integer :: nAtomsInCol
329 – #endif
330 –    integer :: alloc_stat
331 –
332 –
333 –    status = 0
334 – #ifdef IS_MPI
246      nAtomsInRow = getNatomsInRow(plan_atom_row)
247      nAtomsInCol = getNatomsInCol(plan_atom_col)
248   #endif
# Line 342 | Line 253 | contains
253         status = -1
254         return
255      end if
256 +
257      if (allocated(rho)) deallocate(rho)
258      allocate(rho(nlocal),stat=alloc_stat)
259      if (alloc_stat /= 0) then
# Line 356 | Line 268 | contains
268         return
269      end if
270  
359 –    if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
360 –    allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
361 –    if (alloc_stat /= 0) then
362 –       status = -1
363 –       return
364 –    end if
365 –
271   #ifdef IS_MPI
272  
273      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 372 | Line 277 | contains
277         return
278      end if
279  
375 –
280      if (allocated(frho_row)) deallocate(frho_row)
281      allocate(frho_row(nAtomsInRow),stat=alloc_stat)
282      if (alloc_stat /= 0) then
# Line 391 | Line 295 | contains
295         status = -1
296         return
297      end if
394 –    if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
395 –    allocate(d2frhodrhodrho_row(nAtomsInRow),stat=alloc_stat)
396 –    if (alloc_stat /= 0) then
397 –       status = -1
398 –       return
399 –    end if
298  
401 –
299      ! Now do column arrays
300  
301      if (allocated(frho_col)) deallocate(frho_col)
# Line 419 | Line 316 | contains
316         status = -1
317         return
318      end if
422 –    if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
423 –    allocate(d2frhodrhodrho_col(nAtomsInCol),stat=alloc_stat)
424 –    if (alloc_stat /= 0) then
425 –       status = -1
426 –       return
427 –    end if
319  
320   #endif
321  
322 <  end subroutine allocateEAM
322 >  end subroutine init_EAM_FF
323  
324 <  !! C sets rcut to be the largest cutoff of any atype
434 <  !! present in this simulation. Doesn't include all atypes
435 <  !! sim knows about, just those in the simulation.
436 <  subroutine setCutoffEAM(rcut, status)
324 >  subroutine setCutoffEAM(rcut)
325      real(kind=dp) :: rcut
438 –    integer :: status
439 –    status = 0
440 –
326      EAM_rcut = rcut
442 –
327    end subroutine setCutoffEAM
328  
445 –
446 –
329    subroutine clean_EAM()
330  
331      ! clean non-IS_MPI first
# Line 462 | Line 344 | contains
344   #endif
345    end subroutine clean_EAM
346  
465 –
466 –
467 –  subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
468 –    integer, intent(in)          :: eam_n_rho
469 –    integer, intent(in)          :: eam_n_r
470 –    type (EAMType)               :: thisEAMType
471 –    integer, optional   :: stat
472 –    integer             :: alloc_stat
473 –
474 –
475 –
476 –    if (present(stat)) stat = 0
477 –
478 –    allocate(thisEAMType%eam_rvals(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_rhovals(eam_n_rho),stat=alloc_stat)  
484 –    if (alloc_stat /= 0 ) then
485 –       if (present(stat)) stat = -1
486 –       return
487 –    end if
488 –    allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)  
489 –    if (alloc_stat /= 0 ) then
490 –       if (present(stat)) stat = -1
491 –       return
492 –    end if
493 –    allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)        
494 –    if (alloc_stat /= 0 ) then
495 –       if (present(stat)) stat = -1
496 –       return
497 –    end if
498 –    allocate(thisEAMType%eam_rho_r(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_phi_r(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_F_rho_pp(eam_n_rho),stat=alloc_stat)  
509 –    if (alloc_stat /= 0 ) then
510 –       if (present(stat)) stat = -1
511 –       return
512 –    end if
513 –    allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)  
514 –    if (alloc_stat /= 0 ) then
515 –       if (present(stat)) stat = -1
516 –       return
517 –    end if
518 –    allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)  
519 –    if (alloc_stat /= 0 ) then
520 –       if (present(stat)) stat = -1
521 –       return
522 –    end if
523 –    allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
524 –    if (alloc_stat /= 0 ) then
525 –       if (present(stat)) stat = -1
526 –       return
527 –    end if
528 –
529 –
530 –  end subroutine allocate_EAMType
531 –
532 –
347    subroutine deallocate_EAMType(thisEAMType)
348      type (EAMtype), pointer :: thisEAMType
349  
350 <    ! free Arrays in reverse order of allocation...
351 <    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
352 <    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
353 <    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
540 <    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
541 <    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
542 <    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
543 <    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
544 <    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
545 <    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
546 <    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
350 >    call deleteSpline(thisEAMType%F)
351 >    call deleteSpline(thisEAMType%rho)
352 >    call deleteSpline(thisEAMType%phi)
353 >    call deleteSpline(thisEAMType%Z)
354  
355    end subroutine deallocate_EAMType
356  
357    !! Calculates rho_r
358    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
359 <    integer :: atom1,atom2
359 >    integer :: atom1, atom2
360      real(kind = dp), dimension(3) :: d
361      real(kind = dp), intent(inout)               :: r
362      real(kind = dp), intent(inout)               :: rijsq
# Line 557 | Line 364 | contains
364      real(kind = dp) :: rho_i_at_j
365      ! value of electron density rho do to atom j at atom i
366      real(kind = dp) :: rho_j_at_i
560 –
561 –    ! we don't use the derivatives, dummy variables
562 –    real( kind = dp) :: drho,d2rho
367      integer :: eam_err
368      
369 <    integer :: atid1,atid2 ! Global atid    
369 >    integer :: atid1, atid2 ! Global atid    
370      integer :: myid_atom1 ! EAM atid
371      integer :: myid_atom2
372  
569 –
373      ! check to see if we need to be cleaned at the start of a force loop
374  
572 –
573 –
574 –
375   #ifdef IS_MPI
376      Atid1 = Atid_row(Atom1)
377      Atid2 = Atid_col(Atom2)
# Line 585 | Line 385 | contains
385  
386      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
387  
388 <
389 <
590 <       call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
591 <            EAMList%EAMParams(myid_atom1)%eam_rvals, &
592 <            EAMList%EAMParams(myid_atom1)%eam_rho_r, &
593 <            EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
594 <            r, rho_i_at_j,drho,d2rho)
388 >       call lookupUniformSpline(EAMList%EAMParams(myid_atom1)%rho, r, &
389 >            rho_i_at_j)
390  
596 –
597 –
391   #ifdef  IS_MPI
392         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
393   #else
394         rho(atom2) = rho(atom2) + rho_i_at_j
395   #endif
603 –             ! write(*,*) atom1,atom2,r,rho_i_at_j
396      endif
397  
398      if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
607 –       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
608 –            EAMList%EAMParams(myid_atom2)%eam_rvals, &
609 –            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
610 –            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
611 –            r, rho_j_at_i,drho,d2rho)
399  
400 +       call lookupUniformSpline(EAMList%EAMParams(myid_atom2)%rho, r, &
401 +            rho_j_at_i)
402  
614 –
615 –
403   #ifdef  IS_MPI
404         rho_row(atom1) = rho_row(atom1) + rho_j_at_i
405   #else
# Line 620 | Line 407 | contains
407   #endif
408      endif
409  
623 –
624 –
625 –
626 –
627 –
410    end subroutine calc_eam_prepair_rho
411  
412  
631 –
632 –
413    !! Calculate the functional F(rho) for all local atoms
414 <  subroutine calc_eam_preforce_Frho(nlocal,pot)
414 >  subroutine calc_eam_preforce_Frho(nlocal, pot)
415      integer :: nlocal
416      real(kind=dp) :: pot
417 <    integer :: i,j
417 >    integer :: i, j
418      integer :: atom
419 <    real(kind=dp) :: U,U1,U2
419 >    real(kind=dp) :: U,U1
420      integer :: atype1
421 <    integer :: me,atid1
642 <    integer :: n_rho_points
421 >    integer :: me, atid1
422  
644 –
423      cleanme = .true.
424 <    !! Scatter the electron density from  pre-pair calculation back to local atoms
424 >    !! Scatter the electron density from  pre-pair calculation back to
425 >    !! local atoms
426   #ifdef IS_MPI
427      call scatter(rho_row,rho,plan_atom_row,eam_err)
428      if (eam_err /= 0 ) then
# Line 659 | Line 438 | contains
438      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
439   #endif
440  
662 –
663 –
441      !! Calculate F(rho) and derivative
442      do atom = 1, nlocal
443         atid1 = atid(atom)
444         me = eamList%atidToEAMtype(atid1)
668 –       n_rho_points = EAMList%EAMParams(me)%eam_nrho
669 –       !  Check to see that the density is not greater than the larges rho we have calculated
670 –       if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
671 –          call eam_splint(n_rho_points, &
672 –               EAMList%EAMParams(me)%eam_rhovals, &
673 –               EAMList%EAMParams(me)%eam_f_rho, &
674 –               EAMList%EAMParams(me)%eam_f_rho_pp, &
675 –               rho(atom), & ! Actual Rho
676 –               u, u1, u2)
677 –       else
678 –          ! Calculate F(rho with the largest available rho value
679 –          call eam_splint(n_rho_points, &
680 –               EAMList%EAMParams(me)%eam_rhovals, &
681 –               EAMList%EAMParams(me)%eam_f_rho, &
682 –               EAMList%EAMParams(me)%eam_f_rho_pp, &
683 –               EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
684 –               u,u1,u2)
685 –       end if
445  
446 <
446 >       call lookupUniformSpline1d(EAMList%EAMParams(me)%F, rho(atom), &
447 >            u, u1)
448 >      
449         frho(atom) = u
450         dfrhodrho(atom) = u1
690 –       d2frhodrhodrho(atom) = u2
451         pot = pot + u
692 –
452      enddo
453  
695 –
696 –
454   #ifdef IS_MPI
455      !! communicate f(rho) and derivatives back into row and column arrays
456      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 712 | Line 469 | contains
469      if (eam_err /=  0) then
470         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
471      endif
715 –
716 –
717 –
718 –
719 –
720 –    if (nmflag) then
721 –       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row)
722 –       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col)
723 –    endif
472   #endif
473  
474  
475    end subroutine calc_eam_preforce_Frho
476 <
729 <
730 <
731 <
476 >  
477    !! Does EAM pairwise Force calculation.  
478    subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
479         pot, f, do_pot)
# Line 742 | Line 487 | contains
487  
488      logical, intent(in) :: do_pot
489  
490 <    real( kind = dp ) :: drdx,drdy,drdz
491 <    real( kind = dp ) :: d2
492 <    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
493 <    real( kind = dp ) :: rha,drha,d2rha, dpha
749 <    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
490 >    real( kind = dp ) :: drdx, drdy, drdz
491 >    real( kind = dp ) :: phab, pha, dvpdr
492 >    real( kind = dp ) :: rha, drha, dpha
493 >    real( kind = dp ) :: rhb, drhb, dphb
494      real( kind = dp ) :: dudr
495 <    real( kind = dp ) :: rci,rcj
496 <    real( kind = dp ) :: drhoidr,drhojdr
497 <    real( kind = dp ) :: d2rhoidrdr
498 <    real( kind = dp ) :: d2rhojdrdr
755 <    real( kind = dp ) :: Fx,Fy,Fz
756 <    real( kind = dp ) :: r,d2pha,phb,d2phb
495 >    real( kind = dp ) :: rci, rcj
496 >    real( kind = dp ) :: drhoidr, drhojdr
497 >    real( kind = dp ) :: Fx, Fy, Fz
498 >    real( kind = dp ) :: r, phb
499  
500 <    integer :: id1,id2
500 >    integer :: id1, id2
501      integer  :: mytype_atom1
502      integer  :: mytype_atom2
503 <    integer  :: atid1,atid2
762 <    !Local Variables
503 >    integer  :: atid1, atid2
504  
764 –    ! write(*,*) "Frho: ", Frho(atom1)
765 –
505      phab = 0.0E0_DP
506      dvpdr = 0.0E0_DP
768 –    d2vpdrdr = 0.0E0_DP
507  
508      if (rij .lt. EAM_rcut) then
509  
# Line 791 | Line 529 | contains
529         drdz = d(3)/rij
530  
531         if (rij.lt.rci) then
532 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
533 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
534 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
535 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
536 <               rij, rha,drha,d2rha)
537 <          !! Calculate Phi(r) for atom1.
538 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
539 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
540 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
541 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
542 <               rij, pha,dpha,d2pha)
532 >
533 >          ! Calculate rho and drho for atom1
534 >
535 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom1)%rho, &
536 >               rij, rha, drha)
537 >          
538 >          ! Calculate Phi(r) for atom1.
539 >          
540 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom1)%phi, &
541 >               rij, pha, dpha)
542 >
543         endif
544  
545         if (rij.lt.rcj) then      
808 –          ! Calculate rho,drho and d2rho for atom1
809 –          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
810 –               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
811 –               EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
812 –               EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
813 –               rij, rhb,drhb,d2rhb)
546  
547 <          !! Calculate Phi(r) for atom2.
548 <          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
549 <               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
550 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
551 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
552 <               rij, phb,dphb,d2phb)
547 >          ! Calculate rho and drho for atom2
548 >
549 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom2)%rho, &
550 >               rij, rhb, drhb)
551 >
552 >          ! Calculate Phi(r) for atom2.
553 >
554 >          call lookupUniformSpline1d(EAMList%EAMParams(mytype_atom2)%phi, &
555 >               rij, phb, dphb)
556 >
557         endif
558  
559         if (rij.lt.rci) then
560            phab = phab + 0.5E0_DP*(rhb/rha)*pha
561            dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
562                 pha*((drhb/rha) - (rhb*drha/rha/rha)))
827 –          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rhb/rha)*d2pha + &
828 –               2.0E0_DP*dpha*((drhb/rha) - (rhb*drha/rha/rha)) + &
829 –               pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
830 –               (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
563         endif
564  
565         if (rij.lt.rcj) then
566            phab = phab + 0.5E0_DP*(rha/rhb)*phb
567            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
568                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
837 –          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rha/rhb)*d2phb + &
838 –               2.0E0_DP*dphb*((drha/rhb) - (rha*drhb/rhb/rhb)) + &
839 –               phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
840 –               (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
569         endif
570  
571         drhoidr = drha
572         drhojdr = drhb
573  
846 –       d2rhoidrdr = d2rha
847 –       d2rhojdrdr = d2rhb
848 –
849 –
574   #ifdef IS_MPI
575         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
576              + dvpdr
# Line 854 | Line 578 | contains
578   #else
579         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
580              + dvpdr
857 –       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
581   #endif
582  
583         fx = dudr * drdx
# Line 891 | Line 614 | contains
614   #endif
615  
616         vpair = vpair + phab
617 +
618   #ifdef IS_MPI
619         id1 = AtomRowToGlobal(atom1)
620         id2 = AtomColToGlobal(atom2)
# Line 906 | Line 630 | contains
630            fpair(3) = fpair(3) + fz
631  
632         endif
909 –
910 –       if (nmflag) then
911 –
912 –          drhoidr = drha
913 –          drhojdr = drhb
914 –          d2rhoidrdr = d2rha
915 –          d2rhojdrdr = d2rhb
916 –
917 – #ifdef IS_MPI
918 –          d2 = d2vpdrdr + &
919 –               d2rhoidrdr*dfrhodrho_col(atom2) + &
920 –               d2rhojdrdr*dfrhodrho_row(atom1) + &
921 –               drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
922 –               drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
923 –
924 – #else
925 –
926 –          d2 = d2vpdrdr + &
927 –               d2rhoidrdr*dfrhodrho(atom2) + &
928 –               d2rhojdrdr*dfrhodrho(atom1) + &
929 –               drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
930 –               drhojdr*drhojdr*d2frhodrhodrho(atom1)
931 – #endif
932 –       end if
933 –
633      endif
634    end subroutine do_eam_pair
635  
937 –
938 –  subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
939 –
940 –    integer :: atype, nx, j
941 –    real( kind = DP ), dimension(:) :: xa
942 –    real( kind = DP ), dimension(:) :: ya
943 –    real( kind = DP ), dimension(:) :: yppa
944 –    real( kind = DP ) :: x, y
945 –    real( kind = DP ) :: dy, d2y
946 –    real( kind = DP ) :: del, h, a, b, c, d
947 –    integer :: pp_arraySize
948 –
949 –
950 –    ! this spline code assumes that the x points are equally spaced
951 –    ! do not attempt to use this code if they are not.
952 –
953 –
954 –    ! find the closest point with a value below our own:
955 –    j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
956 –
957 –    ! check to make sure we're inside the spline range:
958 –    if ((j.gt.nx).or.(j.lt.1)) then
959 –       write(errMSG,*) "EAM_splint: x is outside bounds of spline: ",x,j
960 –       call handleError(routineName,errMSG)
961 –    endif
962 –    ! check to make sure we haven't screwed up the calculation of j:
963 –    if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
964 –       if (j.ne.nx) then
965 –          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
966 –          call handleError(routineName,errMSG)
967 –       endif
968 –    endif
969 –
970 –    del = xa(j+1) - x
971 –    h = xa(j+1) - xa(j)
972 –
973 –    a = del / h
974 –    b = 1.0E0_DP - a
975 –    c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
976 –    d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
977 –
978 –    y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
979 –
980 –    dy = (ya(j+1)-ya(j))/h &
981 –         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
982 –         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
983 –
984 –
985 –    d2y = a*yppa(j) + b*yppa(j+1)
986 –
987 –
988 –  end subroutine eam_splint
989 –
990 –
991 –  subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
992 –
993 –
994 –    ! yp1 and ypn are the first derivatives of y at the two endpoints
995 –    ! if boundary is 'L' the lower derivative is used
996 –    ! if boundary is 'U' the upper derivative is used
997 –    ! if boundary is 'B' then both derivatives are used
998 –    ! if boundary is anything else, then both derivatives are assumed to be 0
999 –
1000 –    integer :: nx, i, k, max_array_size
1001 –
1002 –    real( kind = DP ), dimension(:)        :: xa
1003 –    real( kind = DP ), dimension(:)        :: ya
1004 –    real( kind = DP ), dimension(:)        :: yppa
1005 –    real( kind = DP ), dimension(size(xa)) :: u
1006 –    real( kind = DP ) :: yp1,ypn,un,qn,sig,p
1007 –    character(len=*) :: boundary
1008 –
1009 –    ! make sure the sizes match
1010 –    if ((nx /= size(xa)) .or. (nx /= size(ya))) then
1011 –       call handleWarning("EAM_SPLINE","Array size mismatch")
1012 –    end if
1013 –
1014 –    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
1015 –         (boundary.eq.'b').or.(boundary.eq.'B')) then
1016 –       yppa(1) = -0.5E0_DP
1017 –       u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
1018 –            ya(1))/(xa(2)-xa(1))-yp1)
1019 –    else
1020 –       yppa(1) = 0.0E0_DP
1021 –       u(1)  = 0.0E0_DP
1022 –    endif
1023 –
1024 –    do i = 2, nx - 1
1025 –       sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1026 –       p = sig * yppa(i-1) + 2.0E0_DP
1027 –       yppa(i) = (sig - 1.0E0_DP) / p
1028 –       u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
1029 –            (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1030 –            (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1031 –    enddo
1032 –
1033 –    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1034 –         (boundary.eq.'b').or.(boundary.eq.'B')) then
1035 –       qn = 0.5E0_DP
1036 –       un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
1037 –            (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
1038 –    else
1039 –       qn = 0.0E0_DP
1040 –       un = 0.0E0_DP
1041 –    endif
1042 –
1043 –    yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1044 –
1045 –    do k = nx-1, 1, -1
1046 –       yppa(k)=yppa(k)*yppa(k+1)+u(k)
1047 –    enddo
1048 –
1049 –  end subroutine eam_spline
1050 –
636   end module eam

Diff Legend

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