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 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2291 by chuckv, Wed Sep 14 20:31:38 2005 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   #ifdef IS_MPI
49    use mpiSimulation
50   #endif
# Line 63 | Line 63 | module eam
63  
64    character(len = 200) :: errMsg
65    character(len=*), parameter :: RoutineName =  "EAM MODULE"
66 < !! Logical that determines if eam arrays should be zeroed
66 >  !! Logical that determines if eam arrays should be zeroed
67    logical :: cleanme = .true.
68    logical :: nmflag  = .false.
69  
70 <
70 >
71    type, private :: EAMtype
72       integer           :: eam_atype      
73       real( kind = DP ) :: eam_dr          
# Line 77 | Line 77 | module eam
77       real( kind = DP ) :: eam_drho      
78       real( kind = DP ) :: eam_rcut    
79       integer           :: eam_atype_map
80 <    
80 >
81       real( kind = DP ), pointer, dimension(:) :: eam_rvals        => null()
82       real( kind = DP ), pointer, dimension(:) :: eam_rhovals      => null()
83       real( kind = DP ), pointer, dimension(:) :: eam_F_rho        => null()
# Line 99 | Line 99 | module eam
99    real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
100  
101  
102 < !! Arrays for MPI storage
102 >  !! Arrays for MPI storage
103   #ifdef IS_MPI
104    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
105    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
# Line 115 | Line 115 | module eam
115    type, private :: EAMTypeList
116       integer           :: n_eam_types = 0
117       integer           :: currentAddition = 0
118 <    
118 >
119       type (EAMtype), pointer  :: EAMParams(:) => null()
120 +     integer, pointer         :: atidToEAMType(:) => null()
121    end type EAMTypeList
122  
123  
# Line 132 | Line 133 | module eam
133    public :: calc_eam_prepair_rho
134    public :: calc_eam_preforce_Frho
135    public :: clean_EAM
136 +  public :: destroyEAMTypes
137 +  public :: getEAMCut
138  
139   contains
140  
141  
142    subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
143         eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
144 <       eam_ident,status)
144 >       c_ident,status)
145      real (kind = dp )                      :: lattice_constant
146      integer                                :: eam_nrho
147      real (kind = dp )                      :: eam_drho
# Line 148 | Line 151 | contains
151      real (kind = dp ), dimension(eam_nr)   :: eam_Z_r
152      real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
153      real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
154 <    integer                                :: eam_ident
154 >    integer                                :: c_ident
155      integer                                :: status
156  
157 <    integer                                :: nAtypes
157 >    integer                                :: nAtypes,nEAMTypes,myATID
158      integer                                :: maxVals
159      integer                                :: alloc_stat
160      integer                                :: current
# Line 162 | Line 165 | contains
165  
166      !! Assume that atypes has already been set and get the total number of types in atypes
167      !! Also assume that every member of atypes is a EAM model.
165 –  
168  
169 +
170      ! check to see if this is the first time into
171      if (.not.associated(EAMList%EAMParams)) then
172 <       call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
173 <       EAMList%n_eam_types = nAtypes
174 <       allocate(EAMList%EAMParams(nAtypes))
172 >       call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
173 >       EAMList%n_eam_types = nEAMtypes
174 >       allocate(EAMList%EAMParams(nEAMTypes))
175 >       nAtypes = getSize(atypes)
176 >       allocate(EAMList%atidToEAMType(nAtypes))
177      end if
178  
179      EAMList%currentAddition = EAMList%currentAddition + 1
180      current = EAMList%currentAddition
176 –    
181  
182 +    myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
183 +    EAMList%atidToEAMType(myATID) = current
184 +
185      call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
186      if (alloc_stat /= 0) then
187         status = -1
188         return
189      end if
190  
191 <    ! this is a possible bug, we assume a correspondence between the vector atypes and
192 <    ! EAMAtypes
186 <      
187 <    EAMList%EAMParams(current)%eam_atype    = eam_ident
191 >  
192 >    EAMList%EAMParams(current)%eam_atype    = c_ident
193      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
194      EAMList%EAMParams(current)%eam_nrho     = eam_nrho
195      EAMList%EAMParams(current)%eam_drho     = eam_drho
# Line 198 | Line 203 | contains
203    end subroutine newEAMtype
204  
205  
206 +  ! kills all eam types entered and sets simulation to uninitalized
207 +  subroutine destroyEAMtypes()
208 +    integer :: i
209 +    type(EAMType), pointer :: tempEAMType=>null()
210 +
211 +    do i = 1, EAMList%n_eam_types
212 +       tempEAMType => eamList%EAMParams(i)
213 +       call deallocate_EAMType(tempEAMType)
214 +    end do
215 +    if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
216 +    eamList%EAMParams => null()
217 +
218 +    eamList%n_eam_types = 0
219 +    eamList%currentAddition = 0
220 +
221 +  end subroutine destroyEAMtypes
222 +
223 +  function getEAMCut(atomID) result(cutValue)
224 +    integer, intent(in) :: atomID
225 +    integer :: eamID
226 +    real(kind=dp) :: cutValue
227 +    
228 +    eamID = EAMList%atidToEAMType(atomID)
229 +    cutValue = EAMList%EAMParams(eamID)%eam_rcut
230 +  end function getEAMCut
231  
232    subroutine init_EAM_FF(status)
233      integer :: status
# Line 214 | Line 244 | contains
244         return
245      end if
246  
217 –
218 –       do i = 1, EAMList%currentAddition
247  
248 < ! Build array of r values
248 >    do i = 1, EAMList%currentAddition
249  
250 <          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
250 >       ! Build array of r values
251  
252 <          ! precompute the pair potential and get it into kcal / mol:
253 <          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
254 <          do j = 2, EAMList%EAMParams(i)%eam_nr
255 <             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
252 >       do j = 1,EAMList%EAMParams(i)%eam_nr
253 >          EAMList%EAMParams(i)%eam_rvals(j) = &
254 >               real(j-1,kind=dp)* &
255 >               EAMList%EAMParams(i)%eam_dr
256         end do
257 <      
257 >       ! Build array of rho values
258 >       do j = 1,EAMList%EAMParams(i)%eam_nrho
259 >          EAMList%EAMParams(i)%eam_rhovals(j) = &
260 >               real(j-1,kind=dp)* &
261 >               EAMList%EAMParams(i)%eam_drho
262 >       end do
263 >       ! convert from eV to kcal / mol:
264 >       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
265  
266 <       do i = 1,  EAMList%currentAddition
267 <          number_r   = EAMList%EAMParams(i)%eam_nr
268 <          number_rho = EAMList%EAMParams(i)%eam_nrho
269 <          
270 <          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')
266 >       ! precompute the pair potential and get it into kcal / mol:
267 >       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
268 >       do j = 2, EAMList%EAMParams(i)%eam_nr
269 >          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
270 >          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
271         enddo
272 +    end do
273  
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
274  
275 < !       EAM_rcut = current_rcut_max
276 < !       EAM_rcut_orig = current_rcut_max
277 < !       do i = 1, EAMList%currentAddition
278 < !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
279 < !       end do
280 <       !! Allocate arrays for force calculation
281 <      
282 <       call allocateEAM(alloc_stat)
283 <       if (alloc_stat /= 0 ) then
284 <          write(*,*) "allocateEAM failed"
285 <          status = -1
286 <          return
287 <       endif
288 <    
275 >    do i = 1,  EAMList%currentAddition
276 >       number_r   = EAMList%EAMParams(i)%eam_nr
277 >       number_rho = EAMList%EAMParams(i)%eam_nrho
278 >
279 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
280 >            EAMList%EAMParams(i)%eam_rho_r, &
281 >            EAMList%EAMParams(i)%eam_rho_r_pp, &
282 >            0.0E0_DP, 0.0E0_DP, 'N')
283 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
284 >            EAMList%EAMParams(i)%eam_Z_r, &
285 >            EAMList%EAMParams(i)%eam_Z_r_pp, &
286 >            0.0E0_DP, 0.0E0_DP, 'N')
287 >       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
288 >            EAMList%EAMParams(i)%eam_F_rho, &
289 >            EAMList%EAMParams(i)%eam_F_rho_pp, &
290 >            0.0E0_DP, 0.0E0_DP, 'N')
291 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
292 >            EAMList%EAMParams(i)%eam_phi_r, &
293 >            EAMList%EAMParams(i)%eam_phi_r_pp, &
294 >            0.0E0_DP, 0.0E0_DP, 'N')
295 >    enddo
296 >
297 >    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
298 >    !! find the smallest rcut for any eam atype
299 >    !       do i = 2, EAMList%currentAddition
300 >    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
301 >    !       end do
302 >
303 >    !       EAM_rcut = current_rcut_max
304 >    !       EAM_rcut_orig = current_rcut_max
305 >    !       do i = 1, EAMList%currentAddition
306 >    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
307 >    !       end do
308 >    !! Allocate arrays for force calculation
309 >
310 >    call allocateEAM(alloc_stat)
311 >    if (alloc_stat /= 0 ) then
312 >       write(*,*) "allocateEAM failed"
313 >       status = -1
314 >       return
315 >    endif
316 >
317    end subroutine init_EAM_FF
318  
319 < !! routine checks to see if array is allocated, deallocates array if allocated
320 < !! and then creates the array to the required size
319 >  !! routine checks to see if array is allocated, deallocates array if allocated
320 >  !! and then creates the array to the required size
321    subroutine allocateEAM(status)
322      integer, intent(out) :: status
323  
# Line 330 | Line 360 | contains
360         status = -1
361         return
362      end if
363 <    
363 >
364   #ifdef IS_MPI
365  
366      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 367 | Line 397 | contains
397      end if
398  
399  
400 < ! Now do column arrays
400 >    ! Now do column arrays
401  
402      if (allocated(frho_col)) deallocate(frho_col)
403      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
# Line 393 | Line 423 | contains
423         status = -1
424         return
425      end if
426 <  
426 >
427   #endif
428  
429    end subroutine allocateEAM
430  
431 < !! C sets rcut to be the largest cutoff of any atype
432 < !! present in this simulation. Doesn't include all atypes
433 < !! sim knows about, just those in the simulation.
431 >  !! C sets rcut to be the largest cutoff of any atype
432 >  !! present in this simulation. Doesn't include all atypes
433 >  !! sim knows about, just those in the simulation.
434    subroutine setCutoffEAM(rcut, status)
435      real(kind=dp) :: rcut
436      integer :: status
# Line 413 | Line 443 | contains
443  
444  
445    subroutine clean_EAM()
446 <  
447 <   ! clean non-IS_MPI first
446 >
447 >    ! clean non-IS_MPI first
448      frho = 0.0_dp
449      rho  = 0.0_dp
450      dfrhodrho = 0.0_dp
451 < ! clean MPI if needed
451 >    ! clean MPI if needed
452   #ifdef IS_MPI
453      frho_row = 0.0_dp
454      frho_col = 0.0_dp
# Line 442 | Line 472 | contains
472  
473  
474      if (present(stat)) stat = 0
475 <    
475 >
476      allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
477      if (alloc_stat /= 0 ) then
478         if (present(stat)) stat = -1
# Line 493 | Line 523 | contains
523         if (present(stat)) stat = -1
524         return
525      end if
496 –      
526  
527 +
528    end subroutine allocate_EAMType
529  
530  
# Line 502 | Line 532 | contains
532      type (EAMtype), pointer :: thisEAMType
533  
534      ! free Arrays in reverse order of allocation...
535 <    deallocate(thisEAMType%eam_phi_r_pp)      
536 <    deallocate(thisEAMType%eam_rho_r_pp)  
537 <    deallocate(thisEAMType%eam_Z_r_pp)  
538 <    deallocate(thisEAMType%eam_F_rho_pp)  
539 <    deallocate(thisEAMType%eam_phi_r)      
540 <    deallocate(thisEAMType%eam_rho_r)      
541 <    deallocate(thisEAMType%eam_Z_r)  
542 <    deallocate(thisEAMType%eam_F_rho)
543 <    deallocate(thisEAMType%eam_rhovals)
544 <    deallocate(thisEAMType%eam_rvals)
545 <  
535 >    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
536 >    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
537 >    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
538 >    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
539 >    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
540 >    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
541 >    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
542 >    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
543 >    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
544 >    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
545 >
546    end subroutine deallocate_EAMType
547  
548 < !! Calculates rho_r
548 >  !! Calculates rho_r
549    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
550      integer :: atom1,atom2
551      real(kind = dp), dimension(3) :: d
# Line 529 | Line 559 | contains
559      ! we don't use the derivatives, dummy variables
560      real( kind = dp) :: drho,d2rho
561      integer :: eam_err
562 <  
563 <    integer :: myid_atom1
564 <    integer :: myid_atom2
562 >    
563 >    integer :: atid1,atid2 ! Global atid    
564 >    integer :: myid_atom1 ! EAM atid
565 >    integer :: myid_atom2
566  
536 – ! check to see if we need to be cleaned at the start of a force loop
537 –      
567  
568 +    ! check to see if we need to be cleaned at the start of a force loop
569  
570  
571 +
572 +
573   #ifdef IS_MPI
574 <    myid_atom1 = atid_Row(atom1)
575 <    myid_atom2 = atid_Col(atom2)
574 >    Atid1 = Atid_row(Atom1)
575 >    Atid2 = Atid_col(Atom2)
576   #else
577 <    myid_atom1 = atid(atom1)
578 <    myid_atom2 = atid(atom2)
577 >    Atid1 = Atid(Atom1)
578 >    Atid2 = Atid(Atom2)
579   #endif
580  
581 +    Myid_atom1 = Eamlist%atidtoeamtype(Atid1)
582 +    Myid_atom2 = Eamlist%atidtoeamtype(Atid2)
583 +
584      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
585  
586  
# Line 557 | Line 592 | contains
592              r, rho_i_at_j,drho,d2rho)
593  
594  
595 <      
595 >
596   #ifdef  IS_MPI
597         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
598   #else
599         rho(atom2) = rho(atom2) + rho_i_at_j
600   #endif
601 < !       write(*,*) atom1,atom2,r,rho_i_at_j
602 <       endif
601 >             ! write(*,*) atom1,atom2,r,rho_i_at_j
602 >    endif
603  
604 <       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
605 <          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
606 <               EAMList%EAMParams(myid_atom2)%eam_rvals, &
607 <               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
608 <               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
609 <               r, rho_j_at_i,drho,d2rho)
604 >    if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
605 >       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
606 >            EAMList%EAMParams(myid_atom2)%eam_rvals, &
607 >            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
608 >            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
609 >            r, rho_j_at_i,drho,d2rho)
610  
611 <    
612 <      
613 <      
611 >
612 >
613 >
614   #ifdef  IS_MPI
615 <          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
615 >       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
616   #else
617 <          rho(atom1) = rho(atom1) + rho_j_at_i
617 >       rho(atom1) = rho(atom1) + rho_j_at_i
618   #endif
619 <       endif
619 >    endif
620  
621  
622  
# Line 601 | Line 636 | contains
636      integer :: atom
637      real(kind=dp) :: U,U1,U2
638      integer :: atype1
639 <    integer :: me
639 >    integer :: me,atid1
640      integer :: n_rho_points
641  
642 <  
642 >
643      cleanme = .true.
644      !! Scatter the electron density from  pre-pair calculation back to local atoms
645   #ifdef IS_MPI
646      call scatter(rho_row,rho,plan_atom_row,eam_err)
647      if (eam_err /= 0 ) then
648 <      write(errMsg,*) " Error scattering rho_row into rho"
649 <      call handleError(RoutineName,errMesg)
650 <   endif      
648 >       write(errMsg,*) " Error scattering rho_row into rho"
649 >       call handleError(RoutineName,errMesg)
650 >    endif
651      call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
652      if (eam_err /= 0 ) then
653 <      write(errMsg,*) " Error scattering rho_col into rho"
654 <      call handleError(RoutineName,errMesg)
655 <   endif
653 >       write(errMsg,*) " Error scattering rho_col into rho"
654 >       call handleError(RoutineName,errMesg)
655 >    endif
656  
657 <      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
657 >    rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
658   #endif
659  
660  
661  
662 < !! Calculate F(rho) and derivative
662 >    !! Calculate F(rho) and derivative
663      do atom = 1, nlocal
664 <       me = atid(atom)
664 >       atid1 = atid(atom)
665 >       me = eamList%atidToEAMtype(atid1)
666         n_rho_points = EAMList%EAMParams(me)%eam_nrho
667         !  Check to see that the density is not greater than the larges rho we have calculated
668         if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
# Line 654 | Line 690 | contains
690  
691      enddo
692  
657 –  
693  
694 +
695   #ifdef IS_MPI
696      !! communicate f(rho) and derivatives back into row and column arrays
697      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 685 | Line 721 | contains
721      endif
722   #endif
723  
724 <  
724 >
725    end subroutine calc_eam_preforce_Frho
726  
727  
# Line 703 | Line 739 | contains
739      real( kind = dp ), intent(inout), dimension(3) :: fpair
740  
741      logical, intent(in) :: do_pot
742 <    
742 >
743      real( kind = dp ) :: drdx,drdy,drdz
744      real( kind = dp ) :: d2
745      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
# Line 720 | Line 756 | contains
756      integer :: id1,id2
757      integer  :: mytype_atom1
758      integer  :: mytype_atom2
759 +    integer  :: atid1,atid2
760 +    !Local Variables
761  
762 < !Local Variables
725 <    
726 <   ! write(*,*) "Frho: ", Frho(atom1)
762 >    ! write(*,*) "Frho: ", Frho(atom1)
763  
764      phab = 0.0E0_DP
765      dvpdr = 0.0E0_DP
# Line 732 | Line 768 | contains
768      if (rij .lt. EAM_rcut) then
769  
770   #ifdef IS_MPI
771 <       mytype_atom1 = atid_row(atom1)
772 <       mytype_atom2 = atid_col(atom2)
771 >       atid1 = atid_row(atom1)
772 >       atid2 = atid_col(atom2)
773   #else
774 <       mytype_atom1 = atid(atom1)
775 <       mytype_atom2 = atid(atom2)
774 >       atid1 = atid(atom1)
775 >       atid2 = atid(atom2)
776   #endif
777 +
778 +       mytype_atom1 = EAMList%atidToEAMType(atid1)
779 +       mytype_atom2 = EAMList%atidTOEAMType(atid2)
780 +
781 +
782         ! get cutoff for atom 1
783         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
784         ! get type specific cutoff for atom 2
785         rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
786 <      
786 >
787         drdx = d(1)/rij
788         drdy = d(2)/rij
789         drdz = d(3)/rij
790 <      
790 >
791         if (rij.lt.rci) then
792            call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
793                 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
# Line 768 | Line 809 | contains
809                 EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
810                 EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
811                 rij, rhb,drhb,d2rhb)
812 <          
812 >
813            !! Calculate Phi(r) for atom2.
814            call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
815                 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
# Line 786 | Line 827 | contains
827                 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
828                 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
829         endif
830 <      
830 >
831         if (rij.lt.rcj) then
832            phab = phab + 0.5E0_DP*(rha/rhb)*phb
833            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 796 | Line 837 | contains
837                 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
838                 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
839         endif
840 <      
840 >
841         drhoidr = drha
842         drhojdr = drhb
843  
# Line 811 | Line 852 | contains
852   #else
853         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
854              + dvpdr
855 <      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
855 >       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
856   #endif
857  
858         fx = dudr * drdx
# Line 828 | Line 869 | contains
869         f_Row(1,atom1) = f_Row(1,atom1) + fx
870         f_Row(2,atom1) = f_Row(2,atom1) + fy
871         f_Row(3,atom1) = f_Row(3,atom1) + fz
872 <      
872 >
873         f_Col(1,atom2) = f_Col(1,atom2) - fx
874         f_Col(2,atom2) = f_Col(2,atom2) - fy
875         f_Col(3,atom2) = f_Col(3,atom2) - fz
# Line 841 | Line 882 | contains
882         f(1,atom1) = f(1,atom1) + fx
883         f(2,atom1) = f(2,atom1) + fy
884         f(3,atom1) = f(3,atom1) + fz
885 <      
885 >
886         f(1,atom2) = f(1,atom2) - fx
887         f(2,atom2) = f(2,atom2) - fy
888         f(3,atom2) = f(3,atom2) - fz
889   #endif
890 <      
890 >
891         vpair = vpair + phab
892   #ifdef IS_MPI
893         id1 = AtomRowToGlobal(atom1)
# Line 855 | Line 896 | contains
896         id1 = atom1
897         id2 = atom2
898   #endif
899 <      
899 >
900         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
901 <          
901 >
902            fpair(1) = fpair(1) + fx
903            fpair(2) = fpair(2) + fy
904            fpair(3) = fpair(3) + fz
905 <          
905 >
906         endif
907 <    
907 >
908         if (nmflag) then
909  
910            drhoidr = drha
# Line 877 | Line 918 | contains
918                 d2rhojdrdr*dfrhodrho_row(atom1) + &
919                 drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
920                 drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
921 <              
921 >
922   #else
923  
924            d2 = d2vpdrdr + &
# Line 887 | Line 928 | contains
928                 drhojdr*drhojdr*d2frhodrhodrho(atom1)
929   #endif
930         end if
931 <      
932 <    endif    
931 >
932 >    endif
933    end subroutine do_eam_pair
934  
935  
# Line 903 | Line 944 | contains
944      real( kind = DP ) :: del, h, a, b, c, d
945      integer :: pp_arraySize
946  
947 <
947 >
948      ! this spline code assumes that the x points are equally spaced
949      ! do not attempt to use this code if they are not.
950 <    
951 <    
950 >
951 >
952      ! find the closest point with a value below our own:
953      j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
954  
# Line 919 | Line 960 | contains
960      ! check to make sure we haven't screwed up the calculation of j:
961      if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
962         if (j.ne.nx) then
963 <        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
964 <       call handleError(routineName,errMSG)
963 >          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
964 >          call handleError(routineName,errMSG)
965         endif
966      endif
967  
968      del = xa(j+1) - x
969      h = xa(j+1) - xa(j)
970 <    
970 >
971      a = del / h
972      b = 1.0E0_DP - a
973      c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
974      d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
975 <    
975 >
976      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 –  
977  
978 +    dy = (ya(j+1)-ya(j))/h &
979 +         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
980 +         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
981 +
982 +
983 +    d2y = a*yppa(j) + b*yppa(j+1)
984 +
985 +
986    end subroutine eam_splint
987  
988  
# Line 953 | Line 994 | contains
994      ! if boundary is 'U' the upper derivative is used
995      ! if boundary is 'B' then both derivatives are used
996      ! if boundary is anything else, then both derivatives are assumed to be 0
997 <    
997 >
998      integer :: nx, i, k, max_array_size
999 <    
999 >
1000      real( kind = DP ), dimension(:)        :: xa
1001      real( kind = DP ), dimension(:)        :: ya
1002      real( kind = DP ), dimension(:)        :: yppa
1003      real( kind = DP ), dimension(size(xa)) :: u
1004      real( kind = DP ) :: yp1,ypn,un,qn,sig,p
1005      character(len=*) :: boundary
1006 <    
1006 >
1007      ! make sure the sizes match
1008      if ((nx /= size(xa)) .or. (nx /= size(ya))) then
1009         call handleWarning("EAM_SPLINE","Array size mismatch")
# Line 977 | Line 1018 | contains
1018         yppa(1) = 0.0E0_DP
1019         u(1)  = 0.0E0_DP
1020      endif
1021 <    
1021 >
1022      do i = 2, nx - 1
1023         sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1024         p = sig * yppa(i-1) + 2.0E0_DP
# Line 986 | Line 1027 | contains
1027              (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1028              (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1029      enddo
1030 <    
1030 >
1031      if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1032           (boundary.eq.'b').or.(boundary.eq.'B')) then
1033         qn = 0.5E0_DP
# Line 998 | Line 1039 | contains
1039      endif
1040  
1041      yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1042 <    
1042 >
1043      do k = nx-1, 1, -1
1044         yppa(k)=yppa(k)*yppa(k+1)+u(k)
1045      enddo
1046  
1047    end subroutine eam_spline
1048  
1008 –
1009 –
1010 –
1049   end module eam
1012 –
1013 –  subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
1014 –       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
1015 –       eam_ident,status)
1016 –       use definitions, ONLY : dp
1017 –       use eam, ONLY : module_newEAMtype => newEAMtype
1018 –    real (kind = dp )                      :: lattice_constant
1019 –    integer                                :: eam_nrho
1020 –    real (kind = dp )                      :: eam_drho
1021 –    integer                                :: eam_nr
1022 –    real (kind = dp )                     :: eam_dr
1023 –    real (kind = dp )                      :: rcut
1024 –    real (kind = dp ), dimension(eam_nr)    :: eam_Z_r
1025 –    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
1026 –    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
1027 –    integer                                :: eam_ident
1028 –    integer                                :: status
1029 –
1030 –    
1031 –    call module_newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
1032 –       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
1033 –       eam_ident,status)
1034 –      
1035 – end subroutine newEAMtype

Diff Legend

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