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 1624 by chuckv, Thu Oct 21 15:25:30 2004 UTC vs.
Revision 2355 by chuckv, Wed Oct 12 18:59:16 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
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 + #define __FORTRAN90
54 + #include "UseTheForce/DarkSide/fInteractionMap.h"
55  
56    INTEGER, PARAMETER :: DP = selected_real_kind(15)
57  
# Line 22 | Line 65 | module eam
65  
66    character(len = 200) :: errMsg
67    character(len=*), parameter :: RoutineName =  "EAM MODULE"
68 < !! Logical that determines if eam arrays should be zeroed
68 >  !! Logical that determines if eam arrays should be zeroed
69    logical :: cleanme = .true.
70    logical :: nmflag  = .false.
71  
72 <
72 >
73    type, private :: EAMtype
74       integer           :: eam_atype      
75       real( kind = DP ) :: eam_dr          
# Line 36 | Line 79 | module eam
79       real( kind = DP ) :: eam_drho      
80       real( kind = DP ) :: eam_rcut    
81       integer           :: eam_atype_map
82 <    
82 >
83       real( kind = DP ), pointer, dimension(:) :: eam_rvals        => null()
84       real( kind = DP ), pointer, dimension(:) :: eam_rhovals      => null()
85       real( kind = DP ), pointer, dimension(:) :: eam_F_rho        => null()
# Line 58 | Line 101 | module eam
101    real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
102  
103  
104 < !! Arrays for MPI storage
104 >  !! Arrays for MPI storage
105   #ifdef IS_MPI
106    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
107    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
# Line 74 | Line 117 | module eam
117    type, private :: EAMTypeList
118       integer           :: n_eam_types = 0
119       integer           :: currentAddition = 0
120 <    
120 >
121       type (EAMtype), pointer  :: EAMParams(:) => null()
122 +     integer, pointer         :: atidToEAMType(:) => null()
123    end type EAMTypeList
124  
125  
# Line 91 | Line 135 | module eam
135    public :: calc_eam_prepair_rho
136    public :: calc_eam_preforce_Frho
137    public :: clean_EAM
138 +  public :: destroyEAMTypes
139 +  public :: getEAMCut
140  
141   contains
142  
143  
144    subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
145         eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
146 <       eam_ident,status)
146 >       c_ident,status)
147      real (kind = dp )                      :: lattice_constant
148      integer                                :: eam_nrho
149      real (kind = dp )                      :: eam_drho
# Line 107 | Line 153 | contains
153      real (kind = dp ), dimension(eam_nr)   :: eam_Z_r
154      real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
155      real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
156 <    integer                                :: eam_ident
156 >    integer                                :: c_ident
157      integer                                :: status
158  
159 <    integer                                :: nAtypes
159 >    integer                                :: nAtypes,nEAMTypes,myATID
160      integer                                :: maxVals
161      integer                                :: alloc_stat
162      integer                                :: current
# Line 121 | Line 167 | contains
167  
168      !! Assume that atypes has already been set and get the total number of types in atypes
169      !! Also assume that every member of atypes is a EAM model.
124 –  
170  
171 +
172      ! check to see if this is the first time into
173      if (.not.associated(EAMList%EAMParams)) then
174 <       call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
175 <       EAMList%n_eam_types = nAtypes
176 <       allocate(EAMList%EAMParams(nAtypes))
174 >       call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
175 >       EAMList%n_eam_types = nEAMtypes
176 >       allocate(EAMList%EAMParams(nEAMTypes))
177 >       nAtypes = getSize(atypes)
178 >       allocate(EAMList%atidToEAMType(nAtypes))
179      end if
180  
181      EAMList%currentAddition = EAMList%currentAddition + 1
182      current = EAMList%currentAddition
135 –    
183  
184 +    myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
185 +    EAMList%atidToEAMType(myATID) = current
186 +
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 <    ! this is a possible bug, we assume a correspondence between the vector atypes and
194 <    ! EAMAtypes
145 <      
146 <    EAMList%EAMParams(current)%eam_atype    = eam_ident
193 >  
194 >    EAMList%EAMParams(current)%eam_atype    = c_ident
195      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
196      EAMList%EAMParams(current)%eam_nrho     = eam_nrho
197      EAMList%EAMParams(current)%eam_drho     = eam_drho
# Line 157 | Line 205 | contains
205    end subroutine newEAMtype
206  
207  
208 +  ! kills all eam types entered and sets simulation to uninitalized
209 +  subroutine destroyEAMtypes()
210 +    integer :: i
211 +    type(EAMType), pointer :: tempEAMType=>null()
212  
213 +    do i = 1, EAMList%n_eam_types
214 +       tempEAMType => eamList%EAMParams(i)
215 +       call deallocate_EAMType(tempEAMType)
216 +    end do
217 +    if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
218 +    eamList%EAMParams => null()
219 +
220 +    eamList%n_eam_types = 0
221 +    eamList%currentAddition = 0
222 +
223 +  end subroutine destroyEAMtypes
224 +
225 +  function getEAMCut(atomID) result(cutValue)
226 +    integer, intent(in) :: atomID
227 +    integer :: eamID
228 +    real(kind=dp) :: cutValue
229 +    
230 +    eamID = EAMList%atidToEAMType(atomID)
231 +    cutValue = EAMList%EAMParams(eamID)%eam_rcut
232 +  end function getEAMCut
233 +
234    subroutine init_EAM_FF(status)
235      integer :: status
236      integer :: i,j
# Line 173 | Line 246 | contains
246         return
247      end if
248  
176 –
177 –       do i = 1, EAMList%currentAddition
249  
250 < ! Build array of r values
250 >    do i = 1, EAMList%currentAddition
251  
252 <          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
252 >       ! Build array of r values
253  
254 <          ! precompute the pair potential and get it into kcal / mol:
255 <          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
256 <          do j = 2, EAMList%EAMParams(i)%eam_nr
257 <             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
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 <      
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 <       do i = 1,  EAMList%currentAddition
269 <          number_r   = EAMList%EAMParams(i)%eam_nr
270 <          number_rho = EAMList%EAMParams(i)%eam_nrho
271 <          
272 <          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')
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  
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
276  
277 < !       EAM_rcut = current_rcut_max
278 < !       EAM_rcut_orig = current_rcut_max
279 < !       do i = 1, EAMList%currentAddition
280 < !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
281 < !       end do
282 <       !! Allocate arrays for force calculation
283 <      
284 <       call allocateEAM(alloc_stat)
285 <       if (alloc_stat /= 0 ) then
286 <          write(*,*) "allocateEAM failed"
287 <          status = -1
288 <          return
289 <       endif
290 <    
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
310 >    !! Allocate arrays for force calculation
311 >
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
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  
# Line 289 | Line 362 | contains
362         status = -1
363         return
364      end if
365 <    
365 >
366   #ifdef IS_MPI
367  
368      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 326 | Line 399 | contains
399      end if
400  
401  
402 < ! Now do column arrays
402 >    ! Now do column arrays
403  
404      if (allocated(frho_col)) deallocate(frho_col)
405      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
# Line 352 | Line 425 | contains
425         status = -1
426         return
427      end if
428 <  
428 >
429   #endif
430  
431    end subroutine allocateEAM
432  
433 < !! 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.
433 >  !! 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)
437      real(kind=dp) :: rcut
438      integer :: status
# Line 372 | Line 445 | contains
445  
446  
447    subroutine clean_EAM()
448 <  
449 <   ! clean non-IS_MPI first
448 >
449 >    ! clean non-IS_MPI first
450      frho = 0.0_dp
451      rho  = 0.0_dp
452      dfrhodrho = 0.0_dp
453 < ! clean MPI if needed
453 >    ! clean MPI if needed
454   #ifdef IS_MPI
455      frho_row = 0.0_dp
456      frho_col = 0.0_dp
# Line 401 | Line 474 | contains
474  
475  
476      if (present(stat)) stat = 0
477 <    
477 >
478      allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
479      if (alloc_stat /= 0 ) then
480         if (present(stat)) stat = -1
# Line 452 | Line 525 | contains
525         if (present(stat)) stat = -1
526         return
527      end if
455 –      
528  
529 +
530    end subroutine allocate_EAMType
531  
532  
# Line 461 | Line 534 | contains
534      type (EAMtype), pointer :: thisEAMType
535  
536      ! free Arrays in reverse order of allocation...
537 <    deallocate(thisEAMType%eam_phi_r_pp)      
538 <    deallocate(thisEAMType%eam_rho_r_pp)  
539 <    deallocate(thisEAMType%eam_Z_r_pp)  
540 <    deallocate(thisEAMType%eam_F_rho_pp)  
541 <    deallocate(thisEAMType%eam_phi_r)      
542 <    deallocate(thisEAMType%eam_rho_r)      
543 <    deallocate(thisEAMType%eam_Z_r)  
544 <    deallocate(thisEAMType%eam_F_rho)
545 <    deallocate(thisEAMType%eam_rhovals)
546 <    deallocate(thisEAMType%eam_rvals)
547 <  
537 >    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
538 >    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
539 >    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)
547 >
548    end subroutine deallocate_EAMType
549  
550 < !! Calculates rho_r
550 >  !! Calculates rho_r
551    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
552      integer :: atom1,atom2
553      real(kind = dp), dimension(3) :: d
# Line 488 | Line 561 | contains
561      ! we don't use the derivatives, dummy variables
562      real( kind = dp) :: drho,d2rho
563      integer :: eam_err
564 <  
565 <    integer :: myid_atom1
566 <    integer :: myid_atom2
564 >    
565 >    integer :: atid1,atid2 ! Global atid    
566 >    integer :: myid_atom1 ! EAM atid
567 >    integer :: myid_atom2
568  
495 – ! check to see if we need to be cleaned at the start of a force loop
496 –      
569  
570 +    ! check to see if we need to be cleaned at the start of a force loop
571  
572  
573 +
574 +
575   #ifdef IS_MPI
576 <    myid_atom1 = atid_Row(atom1)
577 <    myid_atom2 = atid_Col(atom2)
576 >    Atid1 = Atid_row(Atom1)
577 >    Atid2 = Atid_col(Atom2)
578   #else
579 <    myid_atom1 = atid(atom1)
580 <    myid_atom2 = atid(atom2)
579 >    Atid1 = Atid(Atom1)
580 >    Atid2 = Atid(Atom2)
581   #endif
582  
583 +    Myid_atom1 = Eamlist%atidtoeamtype(Atid1)
584 +    Myid_atom2 = Eamlist%atidtoeamtype(Atid2)
585 +
586      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
587  
588  
# Line 516 | Line 594 | contains
594              r, rho_i_at_j,drho,d2rho)
595  
596  
597 <      
597 >
598   #ifdef  IS_MPI
599         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
600   #else
601         rho(atom2) = rho(atom2) + rho_i_at_j
602   #endif
603 < !       write(*,*) atom1,atom2,r,rho_i_at_j
604 <       endif
603 >             ! write(*,*) atom1,atom2,r,rho_i_at_j
604 >    endif
605  
606 <       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)
606 >    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)
612  
613 <    
614 <      
615 <      
613 >
614 >
615 >
616   #ifdef  IS_MPI
617 <          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
617 >       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
618   #else
619 <          rho(atom1) = rho(atom1) + rho_j_at_i
619 >       rho(atom1) = rho(atom1) + rho_j_at_i
620   #endif
621 <       endif
621 >    endif
622  
623  
624  
# Line 560 | Line 638 | contains
638      integer :: atom
639      real(kind=dp) :: U,U1,U2
640      integer :: atype1
641 <    integer :: me
641 >    integer :: me,atid1
642      integer :: n_rho_points
643  
644 <  
644 >
645      cleanme = .true.
646      !! Scatter the electron density from  pre-pair calculation back to local atoms
647   #ifdef IS_MPI
648      call scatter(rho_row,rho,plan_atom_row,eam_err)
649      if (eam_err /= 0 ) then
650 <      write(errMsg,*) " Error scattering rho_row into rho"
651 <      call handleError(RoutineName,errMesg)
652 <   endif      
650 >       write(errMsg,*) " Error scattering rho_row into rho"
651 >       call handleError(RoutineName,errMesg)
652 >    endif
653      call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
654      if (eam_err /= 0 ) then
655 <      write(errMsg,*) " Error scattering rho_col into rho"
656 <      call handleError(RoutineName,errMesg)
657 <   endif
655 >       write(errMsg,*) " Error scattering rho_col into rho"
656 >       call handleError(RoutineName,errMesg)
657 >    endif
658  
659 <      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
659 >    rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
660   #endif
661  
662  
663  
664 < !! Calculate F(rho) and derivative
664 >    !! Calculate F(rho) and derivative
665      do atom = 1, nlocal
666 <       me = atid(atom)
666 >       atid1 = atid(atom)
667 >       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
# Line 613 | Line 692 | contains
692  
693      enddo
694  
616 –  
695  
696 +
697   #ifdef IS_MPI
698      !! communicate f(rho) and derivatives back into row and column arrays
699      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 644 | Line 723 | contains
723      endif
724   #endif
725  
726 <  
726 >
727    end subroutine calc_eam_preforce_Frho
728  
729  
# Line 662 | Line 741 | contains
741      real( kind = dp ), intent(inout), dimension(3) :: fpair
742  
743      logical, intent(in) :: do_pot
744 <    
744 >
745      real( kind = dp ) :: drdx,drdy,drdz
746      real( kind = dp ) :: d2
747      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
# Line 679 | Line 758 | contains
758      integer :: id1,id2
759      integer  :: mytype_atom1
760      integer  :: mytype_atom2
761 +    integer  :: atid1,atid2
762 +    !Local Variables
763  
764 < !Local Variables
684 <    
685 <   ! write(*,*) "Frho: ", Frho(atom1)
764 >    ! write(*,*) "Frho: ", Frho(atom1)
765  
766      phab = 0.0E0_DP
767      dvpdr = 0.0E0_DP
# Line 691 | Line 770 | contains
770      if (rij .lt. EAM_rcut) then
771  
772   #ifdef IS_MPI
773 <       mytype_atom1 = atid_row(atom1)
774 <       mytype_atom2 = atid_col(atom2)
773 >       atid1 = atid_row(atom1)
774 >       atid2 = atid_col(atom2)
775   #else
776 <       mytype_atom1 = atid(atom1)
777 <       mytype_atom2 = atid(atom2)
776 >       atid1 = atid(atom1)
777 >       atid2 = atid(atom2)
778   #endif
779 +
780 +       mytype_atom1 = EAMList%atidToEAMType(atid1)
781 +       mytype_atom2 = EAMList%atidTOEAMType(atid2)
782 +
783 +
784         ! get cutoff for atom 1
785         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
786         ! get type specific cutoff for atom 2
787         rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
788 <      
788 >
789         drdx = d(1)/rij
790         drdy = d(2)/rij
791         drdz = d(3)/rij
792 <      
792 >
793         if (rij.lt.rci) then
794            call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
795                 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
# Line 727 | Line 811 | contains
811                 EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
812                 EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
813                 rij, rhb,drhb,d2rhb)
814 <          
814 >
815            !! Calculate Phi(r) for atom2.
816            call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
817                 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
# Line 745 | Line 829 | contains
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)))
831         endif
832 <      
832 >
833         if (rij.lt.rcj) then
834            phab = phab + 0.5E0_DP*(rha/rhb)*phb
835            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 755 | Line 839 | contains
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)))
841         endif
842 <      
842 >
843         drhoidr = drha
844         drhojdr = drhb
845  
# Line 770 | Line 854 | contains
854   #else
855         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
856              + dvpdr
857 <      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
857 >       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
858   #endif
859  
860         fx = dudr * drdx
# Line 780 | Line 864 | contains
864  
865   #ifdef IS_MPI
866         if (do_pot) then
867 <          pot_Row(atom1) = pot_Row(atom1) + phab*0.5
868 <          pot_Col(atom2) = pot_Col(atom2) + phab*0.5
867 >          pot_Row(EAM_POT,atom1) = pot_Row(EAM_POT,atom1) + phab*0.5
868 >          pot_Col(EAM_POT,atom2) = pot_Col(EAM_POT,atom2) + phab*0.5
869         end if
870  
871         f_Row(1,atom1) = f_Row(1,atom1) + fx
872         f_Row(2,atom1) = f_Row(2,atom1) + fy
873         f_Row(3,atom1) = f_Row(3,atom1) + fz
874 <      
874 >
875         f_Col(1,atom2) = f_Col(1,atom2) - fx
876         f_Col(2,atom2) = f_Col(2,atom2) - fy
877         f_Col(3,atom2) = f_Col(3,atom2) - fz
# Line 800 | Line 884 | contains
884         f(1,atom1) = f(1,atom1) + fx
885         f(2,atom1) = f(2,atom1) + fy
886         f(3,atom1) = f(3,atom1) + fz
887 <      
887 >
888         f(1,atom2) = f(1,atom2) - fx
889         f(2,atom2) = f(2,atom2) - fy
890         f(3,atom2) = f(3,atom2) - fz
891   #endif
892 <      
892 >
893         vpair = vpair + phab
894   #ifdef IS_MPI
895         id1 = AtomRowToGlobal(atom1)
# Line 814 | Line 898 | contains
898         id1 = atom1
899         id2 = atom2
900   #endif
901 <      
901 >
902         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
903 <          
903 >
904            fpair(1) = fpair(1) + fx
905            fpair(2) = fpair(2) + fy
906            fpair(3) = fpair(3) + fz
907 <          
907 >
908         endif
909 <    
909 >
910         if (nmflag) then
911  
912            drhoidr = drha
# Line 836 | Line 920 | contains
920                 d2rhojdrdr*dfrhodrho_row(atom1) + &
921                 drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
922                 drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
923 <              
923 >
924   #else
925  
926            d2 = d2vpdrdr + &
# Line 846 | Line 930 | contains
930                 drhojdr*drhojdr*d2frhodrhodrho(atom1)
931   #endif
932         end if
933 <      
934 <    endif    
933 >
934 >    endif
935    end subroutine do_eam_pair
936  
937  
# Line 862 | Line 946 | contains
946      real( kind = DP ) :: del, h, a, b, c, d
947      integer :: pp_arraySize
948  
949 <
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 <    
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  
# Line 878 | Line 962 | contains
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)
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 <    
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 <    
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
899 <  
900 <  
901 <       d2y = a*yppa(j) + b*yppa(j+1)
902 <  
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  
# Line 912 | Line 996 | contains
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 <    
999 >
1000      integer :: nx, i, k, max_array_size
1001 <    
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 <    
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")
# Line 936 | Line 1020 | contains
1020         yppa(1) = 0.0E0_DP
1021         u(1)  = 0.0E0_DP
1022      endif
1023 <    
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
# Line 945 | Line 1029 | contains
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 <    
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
# Line 957 | Line 1041 | contains
1041      endif
1042  
1043      yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1044 <    
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  
967 –
968 –
969 –
1051   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