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 1608 by gezelter, Wed Oct 20 04:02:48 2004 UTC vs.
Revision 2278 by chrisfen, Fri Aug 26 22:39:25 2005 UTC

# Line 1 | Line 1
1 + !!
2 + !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 + !!
4 + !! The University of Notre Dame grants you ("Licensee") a
5 + !! non-exclusive, royalty free, license to use, modify and
6 + !! redistribute this software in source and binary code form, provided
7 + !! that the following conditions are met:
8 + !!
9 + !! 1. Acknowledgement of the program authors must be made in any
10 + !!    publication of scientific results based in part on use of the
11 + !!    program.  An acceptable form of acknowledgement is citation of
12 + !!    the article in which the program was described (Matthew
13 + !!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + !!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + !!    Parallel Simulation Engine for Molecular Dynamics,"
16 + !!    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + !!
18 + !! 2. Redistributions of source code must retain the above copyright
19 + !!    notice, this list of conditions and the following disclaimer.
20 + !!
21 + !! 3. Redistributions in binary form must reproduce the above copyright
22 + !!    notice, this list of conditions and the following disclaimer in the
23 + !!    documentation and/or other materials provided with the
24 + !!    distribution.
25 + !!
26 + !! This software is provided "AS IS," without a warranty of any
27 + !! kind. All express or implied conditions, representations and
28 + !! warranties, including any implied warranty of merchantability,
29 + !! fitness for a particular purpose or non-infringement, are hereby
30 + !! excluded.  The University of Notre Dame and its licensors shall not
31 + !! be liable for any damages suffered by licensee as a result of
32 + !! using, modifying or distributing the software or its
33 + !! derivatives. In no event will the University of Notre Dame or its
34 + !! licensors be liable for any lost revenue, profit or data, or for
35 + !! direct, indirect, special, consequential, incidental or punitive
36 + !! damages, however caused and regardless of the theory of liability,
37 + !! arising out of the use of or inability to use software, even if the
38 + !! University of Notre Dame has been advised of the possibility of
39 + !! such damages.
40 + !!
41 +
42   module eam
2 –  use definitions, ONLY : DP,default_error
43    use simulation
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
51    implicit none
52    PRIVATE
53  
54 +  INTEGER, PARAMETER :: DP = selected_real_kind(15)
55  
56    logical, save :: EAM_FF_initialized = .false.
57    integer, save :: EAM_Mixing_Policy
# Line 22 | 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 36 | 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 58 | 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 74 | 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 91 | 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 107 | 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 121 | 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.
124 –  
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
135 –    
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
145 <      
146 <    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 157 | 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 +
231 +  end function getEAMCut
232 +
233    subroutine init_EAM_FF(status)
234      integer :: status
235      integer :: i,j
# Line 173 | Line 245 | contains
245         return
246      end if
247  
176 –
177 –       do i = 1, EAMList%currentAddition
248  
249 < ! Build array of r values
249 >    do i = 1, EAMList%currentAddition
250  
251 <          do j = 1,EAMList%EAMParams(i)%eam_nr
182 <             EAMList%EAMParams(i)%eam_rvals(j) = &
183 <                  real(j-1,kind=dp)* &
184 <                  EAMList%EAMParams(i)%eam_dr
185 <              end do
186 < ! Build array of rho values
187 <          do j = 1,EAMList%EAMParams(i)%eam_nrho
188 <             EAMList%EAMParams(i)%eam_rhovals(j) = &
189 <                  real(j-1,kind=dp)* &
190 <                  EAMList%EAMParams(i)%eam_drho
191 <          end do
192 <          ! convert from eV to kcal / mol:
193 <          EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
251 >       ! Build array of r values
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)
199 <             EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
200 <          enddo
253 >       do j = 1,EAMList%EAMParams(i)%eam_nr
254 >          EAMList%EAMParams(i)%eam_rvals(j) = &
255 >               real(j-1,kind=dp)* &
256 >               EAMList%EAMParams(i)%eam_dr
257         end do
258 <      
258 >       ! Build array of rho values
259 >       do j = 1,EAMList%EAMParams(i)%eam_nrho
260 >          EAMList%EAMParams(i)%eam_rhovals(j) = &
261 >               real(j-1,kind=dp)* &
262 >               EAMList%EAMParams(i)%eam_drho
263 >       end do
264 >       ! convert from eV to kcal / mol:
265 >       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
266  
267 <       do i = 1,  EAMList%currentAddition
268 <          number_r   = EAMList%EAMParams(i)%eam_nr
269 <          number_rho = EAMList%EAMParams(i)%eam_nrho
270 <          
271 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
209 <               EAMList%EAMParams(i)%eam_rho_r, &
210 <               EAMList%EAMParams(i)%eam_rho_r_pp, &
211 <               0.0E0_DP, 0.0E0_DP, 'N')
212 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
213 <               EAMList%EAMParams(i)%eam_Z_r, &
214 <               EAMList%EAMParams(i)%eam_Z_r_pp, &
215 <               0.0E0_DP, 0.0E0_DP, 'N')
216 <          call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
217 <               EAMList%EAMParams(i)%eam_F_rho, &
218 <               EAMList%EAMParams(i)%eam_F_rho_pp, &
219 <               0.0E0_DP, 0.0E0_DP, 'N')
220 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
221 <               EAMList%EAMParams(i)%eam_phi_r, &
222 <               EAMList%EAMParams(i)%eam_phi_r_pp, &
223 <               0.0E0_DP, 0.0E0_DP, 'N')
267 >       ! precompute the pair potential and get it into kcal / mol:
268 >       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
269 >       do j = 2, EAMList%EAMParams(i)%eam_nr
270 >          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
271 >          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
272         enddo
273 +    end do
274  
226 – !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
227 –       !! find the smallest rcut for any eam atype
228 – !       do i = 2, EAMList%currentAddition
229 – !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
230 – !       end do
275  
276 < !       EAM_rcut = current_rcut_max
277 < !       EAM_rcut_orig = current_rcut_max
278 < !       do i = 1, EAMList%currentAddition
279 < !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
280 < !       end do
281 <       !! Allocate arrays for force calculation
282 <      
283 <       call allocateEAM(alloc_stat)
284 <       if (alloc_stat /= 0 ) then
285 <          write(*,*) "allocateEAM failed"
286 <          status = -1
287 <          return
288 <       endif
289 <    
276 >    do i = 1,  EAMList%currentAddition
277 >       number_r   = EAMList%EAMParams(i)%eam_nr
278 >       number_rho = EAMList%EAMParams(i)%eam_nrho
279 >
280 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
281 >            EAMList%EAMParams(i)%eam_rho_r, &
282 >            EAMList%EAMParams(i)%eam_rho_r_pp, &
283 >            0.0E0_DP, 0.0E0_DP, 'N')
284 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
285 >            EAMList%EAMParams(i)%eam_Z_r, &
286 >            EAMList%EAMParams(i)%eam_Z_r_pp, &
287 >            0.0E0_DP, 0.0E0_DP, 'N')
288 >       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
289 >            EAMList%EAMParams(i)%eam_F_rho, &
290 >            EAMList%EAMParams(i)%eam_F_rho_pp, &
291 >            0.0E0_DP, 0.0E0_DP, 'N')
292 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
293 >            EAMList%EAMParams(i)%eam_phi_r, &
294 >            EAMList%EAMParams(i)%eam_phi_r_pp, &
295 >            0.0E0_DP, 0.0E0_DP, 'N')
296 >    enddo
297 >
298 >    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
299 >    !! find the smallest rcut for any eam atype
300 >    !       do i = 2, EAMList%currentAddition
301 >    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
302 >    !       end do
303 >
304 >    !       EAM_rcut = current_rcut_max
305 >    !       EAM_rcut_orig = current_rcut_max
306 >    !       do i = 1, EAMList%currentAddition
307 >    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
308 >    !       end do
309 >    !! Allocate arrays for force calculation
310 >
311 >    call allocateEAM(alloc_stat)
312 >    if (alloc_stat /= 0 ) then
313 >       write(*,*) "allocateEAM failed"
314 >       status = -1
315 >       return
316 >    endif
317 >
318    end subroutine init_EAM_FF
319  
320 < !! routine checks to see if array is allocated, deallocates array if allocated
321 < !! and then creates the array to the required size
320 >  !! routine checks to see if array is allocated, deallocates array if allocated
321 >  !! and then creates the array to the required size
322    subroutine allocateEAM(status)
323      integer, intent(out) :: status
324  
# Line 289 | Line 361 | contains
361         status = -1
362         return
363      end if
364 <    
364 >
365   #ifdef IS_MPI
366  
367      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 326 | Line 398 | contains
398      end if
399  
400  
401 < ! Now do column arrays
401 >    ! Now do column arrays
402  
403      if (allocated(frho_col)) deallocate(frho_col)
404      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
# Line 352 | Line 424 | contains
424         status = -1
425         return
426      end if
427 <  
427 >
428   #endif
429  
430    end subroutine allocateEAM
431  
432 < !! C sets rcut to be the largest cutoff of any atype
433 < !! present in this simulation. Doesn't include all atypes
434 < !! sim knows about, just those in the simulation.
432 >  !! C sets rcut to be the largest cutoff of any atype
433 >  !! present in this simulation. Doesn't include all atypes
434 >  !! sim knows about, just those in the simulation.
435    subroutine setCutoffEAM(rcut, status)
436      real(kind=dp) :: rcut
437      integer :: status
# Line 372 | Line 444 | contains
444  
445  
446    subroutine clean_EAM()
447 <  
448 <   ! clean non-IS_MPI first
447 >
448 >    ! clean non-IS_MPI first
449      frho = 0.0_dp
450      rho  = 0.0_dp
451      dfrhodrho = 0.0_dp
452 < ! clean MPI if needed
452 >    ! clean MPI if needed
453   #ifdef IS_MPI
454      frho_row = 0.0_dp
455      frho_col = 0.0_dp
# Line 401 | Line 473 | contains
473  
474  
475      if (present(stat)) stat = 0
476 <    
476 >
477      allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
478      if (alloc_stat /= 0 ) then
479         if (present(stat)) stat = -1
# Line 452 | Line 524 | contains
524         if (present(stat)) stat = -1
525         return
526      end if
455 –      
527  
528 +
529    end subroutine allocate_EAMType
530  
531  
# Line 461 | Line 533 | contains
533      type (EAMtype), pointer :: thisEAMType
534  
535      ! free Arrays in reverse order of allocation...
536 <    deallocate(thisEAMType%eam_phi_r_pp)      
537 <    deallocate(thisEAMType%eam_rho_r_pp)  
538 <    deallocate(thisEAMType%eam_Z_r_pp)  
539 <    deallocate(thisEAMType%eam_F_rho_pp)  
540 <    deallocate(thisEAMType%eam_phi_r)      
541 <    deallocate(thisEAMType%eam_rho_r)      
542 <    deallocate(thisEAMType%eam_Z_r)  
543 <    deallocate(thisEAMType%eam_F_rho)
544 <    deallocate(thisEAMType%eam_rhovals)
545 <    deallocate(thisEAMType%eam_rvals)
546 <  
536 >    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
537 >    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
538 >    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
539 >    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
540 >    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
541 >    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
542 >    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
543 >    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
544 >    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
545 >    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
546 >
547    end subroutine deallocate_EAMType
548  
549 < !! Calculates rho_r
549 >  !! Calculates rho_r
550    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
551      integer :: atom1,atom2
552      real(kind = dp), dimension(3) :: d
# Line 488 | Line 560 | contains
560      ! we don't use the derivatives, dummy variables
561      real( kind = dp) :: drho,d2rho
562      integer :: eam_err
563 <  
563 >
564      integer :: myid_atom1
565      integer :: myid_atom2
566  
567 < ! check to see if we need to be cleaned at the start of a force loop
496 <      
567 >    ! check to see if we need to be cleaned at the start of a force loop
568  
569  
570  
571 +
572   #ifdef IS_MPI
573      myid_atom1 = atid_Row(atom1)
574      myid_atom2 = atid_Col(atom2)
# Line 516 | Line 588 | contains
588              r, rho_i_at_j,drho,d2rho)
589  
590  
591 <      
591 >
592   #ifdef  IS_MPI
593         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
594   #else
595         rho(atom2) = rho(atom2) + rho_i_at_j
596   #endif
597 < !       write(*,*) atom1,atom2,r,rho_i_at_j
598 <       endif
597 >       !       write(*,*) atom1,atom2,r,rho_i_at_j
598 >    endif
599  
600 <       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
601 <          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
602 <               EAMList%EAMParams(myid_atom2)%eam_rvals, &
603 <               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
604 <               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
605 <               r, rho_j_at_i,drho,d2rho)
600 >    if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
601 >       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
602 >            EAMList%EAMParams(myid_atom2)%eam_rvals, &
603 >            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
604 >            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
605 >            r, rho_j_at_i,drho,d2rho)
606  
607 <    
608 <      
609 <      
607 >
608 >
609 >
610   #ifdef  IS_MPI
611 <          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
611 >       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
612   #else
613 <          rho(atom1) = rho(atom1) + rho_j_at_i
613 >       rho(atom1) = rho(atom1) + rho_j_at_i
614   #endif
615 <       endif
615 >    endif
616  
617  
618  
# Line 563 | Line 635 | contains
635      integer :: me
636      integer :: n_rho_points
637  
638 <  
638 >
639      cleanme = .true.
640      !! Scatter the electron density from  pre-pair calculation back to local atoms
641   #ifdef IS_MPI
642      call scatter(rho_row,rho,plan_atom_row,eam_err)
643      if (eam_err /= 0 ) then
644 <      write(errMsg,*) " Error scattering rho_row into rho"
645 <      call handleError(RoutineName,errMesg)
646 <   endif      
644 >       write(errMsg,*) " Error scattering rho_row into rho"
645 >       call handleError(RoutineName,errMesg)
646 >    endif
647      call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
648      if (eam_err /= 0 ) then
649 <      write(errMsg,*) " Error scattering rho_col into rho"
650 <      call handleError(RoutineName,errMesg)
651 <   endif
649 >       write(errMsg,*) " Error scattering rho_col into rho"
650 >       call handleError(RoutineName,errMesg)
651 >    endif
652  
653 <      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
653 >    rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
654   #endif
655  
656  
657  
658 < !! Calculate F(rho) and derivative
658 >    !! Calculate F(rho) and derivative
659      do atom = 1, nlocal
660         me = atid(atom)
661         n_rho_points = EAMList%EAMParams(me)%eam_nrho
# Line 613 | Line 685 | contains
685  
686      enddo
687  
616 –  
688  
689 +
690   #ifdef IS_MPI
691      !! communicate f(rho) and derivatives back into row and column arrays
692      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 644 | Line 716 | contains
716      endif
717   #endif
718  
719 <  
719 >
720    end subroutine calc_eam_preforce_Frho
721  
722  
# Line 662 | Line 734 | contains
734      real( kind = dp ), intent(inout), dimension(3) :: fpair
735  
736      logical, intent(in) :: do_pot
737 <    
737 >
738      real( kind = dp ) :: drdx,drdy,drdz
739      real( kind = dp ) :: d2
740      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
# Line 680 | Line 752 | contains
752      integer  :: mytype_atom1
753      integer  :: mytype_atom2
754  
755 < !Local Variables
684 <    
685 <   ! write(*,*) "Frho: ", Frho(atom1)
755 >    !Local Variables
756  
757 +    ! write(*,*) "Frho: ", Frho(atom1)
758 +
759      phab = 0.0E0_DP
760      dvpdr = 0.0E0_DP
761      d2vpdrdr = 0.0E0_DP
# Line 701 | Line 773 | contains
773         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
774         ! get type specific cutoff for atom 2
775         rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
776 <      
776 >
777         drdx = d(1)/rij
778         drdy = d(2)/rij
779         drdz = d(3)/rij
780 <      
780 >
781         if (rij.lt.rci) then
782            call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
783                 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
# Line 727 | Line 799 | contains
799                 EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
800                 EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
801                 rij, rhb,drhb,d2rhb)
802 <          
802 >
803            !! Calculate Phi(r) for atom2.
804            call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
805                 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
# Line 745 | Line 817 | contains
817                 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
818                 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
819         endif
820 <      
820 >
821         if (rij.lt.rcj) then
822            phab = phab + 0.5E0_DP*(rha/rhb)*phb
823            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 755 | Line 827 | contains
827                 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
828                 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
829         endif
830 <      
830 >
831         drhoidr = drha
832         drhojdr = drhb
833  
# Line 770 | Line 842 | contains
842   #else
843         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
844              + dvpdr
845 <      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
845 >       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
846   #endif
847  
848         fx = dudr * drdx
# Line 787 | Line 859 | contains
859         f_Row(1,atom1) = f_Row(1,atom1) + fx
860         f_Row(2,atom1) = f_Row(2,atom1) + fy
861         f_Row(3,atom1) = f_Row(3,atom1) + fz
862 <      
862 >
863         f_Col(1,atom2) = f_Col(1,atom2) - fx
864         f_Col(2,atom2) = f_Col(2,atom2) - fy
865         f_Col(3,atom2) = f_Col(3,atom2) - fz
# Line 800 | Line 872 | contains
872         f(1,atom1) = f(1,atom1) + fx
873         f(2,atom1) = f(2,atom1) + fy
874         f(3,atom1) = f(3,atom1) + fz
875 <      
875 >
876         f(1,atom2) = f(1,atom2) - fx
877         f(2,atom2) = f(2,atom2) - fy
878         f(3,atom2) = f(3,atom2) - fz
879   #endif
880 <      
880 >
881         vpair = vpair + phab
882   #ifdef IS_MPI
883         id1 = AtomRowToGlobal(atom1)
# Line 814 | Line 886 | contains
886         id1 = atom1
887         id2 = atom2
888   #endif
889 <      
889 >
890         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
891 <          
891 >
892            fpair(1) = fpair(1) + fx
893            fpair(2) = fpair(2) + fy
894            fpair(3) = fpair(3) + fz
895 <          
895 >
896         endif
897 <    
897 >
898         if (nmflag) then
899  
900            drhoidr = drha
# Line 836 | Line 908 | contains
908                 d2rhojdrdr*dfrhodrho_row(atom1) + &
909                 drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
910                 drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
911 <              
911 >
912   #else
913  
914            d2 = d2vpdrdr + &
# Line 846 | Line 918 | contains
918                 drhojdr*drhojdr*d2frhodrhodrho(atom1)
919   #endif
920         end if
921 <      
922 <    endif    
921 >
922 >    endif
923    end subroutine do_eam_pair
924  
925  
# Line 862 | Line 934 | contains
934      real( kind = DP ) :: del, h, a, b, c, d
935      integer :: pp_arraySize
936  
937 <
937 >
938      ! this spline code assumes that the x points are equally spaced
939      ! do not attempt to use this code if they are not.
940 <    
941 <    
940 >
941 >
942      ! find the closest point with a value below our own:
943      j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
944  
# Line 878 | Line 950 | contains
950      ! check to make sure we haven't screwed up the calculation of j:
951      if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
952         if (j.ne.nx) then
953 <        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
954 <       call handleError(routineName,errMSG)
953 >          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
954 >          call handleError(routineName,errMSG)
955         endif
956      endif
957  
958      del = xa(j+1) - x
959      h = xa(j+1) - xa(j)
960 <    
960 >
961      a = del / h
962      b = 1.0E0_DP - a
963      c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
964      d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
965 <    
965 >
966      y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
895 –  
896 –       dy = (ya(j+1)-ya(j))/h &
897 –            - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
898 –            + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
899 –  
900 –  
901 –       d2y = a*yppa(j) + b*yppa(j+1)
902 –  
967  
968 +    dy = (ya(j+1)-ya(j))/h &
969 +         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
970 +         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
971 +
972 +
973 +    d2y = a*yppa(j) + b*yppa(j+1)
974 +
975 +
976    end subroutine eam_splint
977  
978  
# Line 912 | Line 984 | contains
984      ! if boundary is 'U' the upper derivative is used
985      ! if boundary is 'B' then both derivatives are used
986      ! if boundary is anything else, then both derivatives are assumed to be 0
987 <    
987 >
988      integer :: nx, i, k, max_array_size
989 <    
989 >
990      real( kind = DP ), dimension(:)        :: xa
991      real( kind = DP ), dimension(:)        :: ya
992      real( kind = DP ), dimension(:)        :: yppa
993      real( kind = DP ), dimension(size(xa)) :: u
994      real( kind = DP ) :: yp1,ypn,un,qn,sig,p
995      character(len=*) :: boundary
996 <    
996 >
997      ! make sure the sizes match
998      if ((nx /= size(xa)) .or. (nx /= size(ya))) then
999         call handleWarning("EAM_SPLINE","Array size mismatch")
# Line 936 | Line 1008 | contains
1008         yppa(1) = 0.0E0_DP
1009         u(1)  = 0.0E0_DP
1010      endif
1011 <    
1011 >
1012      do i = 2, nx - 1
1013         sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1014         p = sig * yppa(i-1) + 2.0E0_DP
# Line 945 | Line 1017 | contains
1017              (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1018              (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1019      enddo
1020 <    
1020 >
1021      if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1022           (boundary.eq.'b').or.(boundary.eq.'B')) then
1023         qn = 0.5E0_DP
# Line 957 | Line 1029 | contains
1029      endif
1030  
1031      yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1032 <    
1032 >
1033      do k = nx-1, 1, -1
1034         yppa(k)=yppa(k)*yppa(k+1)+u(k)
1035      enddo
1036  
1037    end subroutine eam_spline
1038  
967 –
968 –
969 –
1039   end module eam
971 –
972 –  subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
973 –       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
974 –       eam_ident,status)
975 –       use definitions, ONLY : dp
976 –       use eam, ONLY : module_newEAMtype => newEAMtype
977 –    real (kind = dp )                      :: lattice_constant
978 –    integer                                :: eam_nrho
979 –    real (kind = dp )                      :: eam_drho
980 –    integer                                :: eam_nr
981 –    real (kind = dp )                     :: eam_dr
982 –    real (kind = dp )                      :: rcut
983 –    real (kind = dp ), dimension(eam_nr)    :: eam_Z_r
984 –    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
985 –    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
986 –    integer                                :: eam_ident
987 –    integer                                :: status
988 –
989 –    
990 –    call module_newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
991 –       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
992 –       eam_ident,status)
993 –      
994 – end subroutine newEAMtype

Diff Legend

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