ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/UseTheForce/DarkSide/eam.F90
(Generate patch)

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/eam.F90 (file contents):
Revision 1608 by gezelter, Wed Oct 20 04:02:48 2004 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 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
# Line 11 | Line 51 | module eam
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    end type EAMTypeList
121  
# Line 91 | Line 132 | module eam
132    public :: calc_eam_prepair_rho
133    public :: calc_eam_preforce_Frho
134    public :: clean_EAM
135 +  public :: destroyEAMTypes
136  
137   contains
138  
# Line 121 | Line 163 | contains
163  
164      !! Assume that atypes has already been set and get the total number of types in atypes
165      !! Also assume that every member of atypes is a EAM model.
124 –  
166  
167 +
168      ! check to see if this is the first time into
169      if (.not.associated(EAMList%EAMParams)) then
170         call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
# Line 132 | Line 174 | contains
174  
175      EAMList%currentAddition = EAMList%currentAddition + 1
176      current = EAMList%currentAddition
135 –    
177  
178 +
179      call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
180      if (alloc_stat /= 0) then
181         status = -1
# Line 142 | Line 184 | contains
184  
185      ! this is a possible bug, we assume a correspondence between the vector atypes and
186      ! EAMAtypes
187 <      
187 >
188      EAMList%EAMParams(current)%eam_atype    = eam_ident
189      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
190      EAMList%EAMParams(current)%eam_nrho     = eam_nrho
# Line 157 | Line 199 | contains
199    end subroutine newEAMtype
200  
201  
202 +  ! kills all eam types entered and sets simulation to uninitalized
203 +  subroutine destroyEAMtypes()
204 +    integer :: i
205 +    type(EAMType), pointer :: tempEAMType=>null()
206  
207 +    do i = 1, EAMList%n_eam_types
208 +       tempEAMType => eamList%EAMParams(i)
209 +       call deallocate_EAMType(tempEAMType)
210 +    end do
211 +    if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
212 +    eamList%EAMParams => null()
213 +
214 +    eamList%n_eam_types = 0
215 +    eamList%currentAddition = 0
216 +
217 +  end subroutine destroyEAMtypes
218 +
219    subroutine init_EAM_FF(status)
220      integer :: status
221      integer :: i,j
# Line 173 | Line 231 | contains
231         return
232      end if
233  
176 –
177 –       do i = 1, EAMList%currentAddition
234  
235 < ! Build array of r values
235 >    do i = 1, EAMList%currentAddition
236  
237 <          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
237 >       ! Build array of r values
238  
239 <          ! precompute the pair potential and get it into kcal / mol:
240 <          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
241 <          do j = 2, EAMList%EAMParams(i)%eam_nr
242 <             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
239 >       do j = 1,EAMList%EAMParams(i)%eam_nr
240 >          EAMList%EAMParams(i)%eam_rvals(j) = &
241 >               real(j-1,kind=dp)* &
242 >               EAMList%EAMParams(i)%eam_dr
243         end do
244 <      
244 >       ! Build array of rho values
245 >       do j = 1,EAMList%EAMParams(i)%eam_nrho
246 >          EAMList%EAMParams(i)%eam_rhovals(j) = &
247 >               real(j-1,kind=dp)* &
248 >               EAMList%EAMParams(i)%eam_drho
249 >       end do
250 >       ! convert from eV to kcal / mol:
251 >       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
252  
253 <       do i = 1,  EAMList%currentAddition
254 <          number_r   = EAMList%EAMParams(i)%eam_nr
255 <          number_rho = EAMList%EAMParams(i)%eam_nrho
256 <          
257 <          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')
253 >       ! precompute the pair potential and get it into kcal / mol:
254 >       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
255 >       do j = 2, EAMList%EAMParams(i)%eam_nr
256 >          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
257 >          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
258         enddo
259 +    end do
260  
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
261  
262 < !       EAM_rcut = current_rcut_max
263 < !       EAM_rcut_orig = current_rcut_max
264 < !       do i = 1, EAMList%currentAddition
265 < !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
266 < !       end do
267 <       !! Allocate arrays for force calculation
268 <      
269 <       call allocateEAM(alloc_stat)
270 <       if (alloc_stat /= 0 ) then
271 <          write(*,*) "allocateEAM failed"
272 <          status = -1
273 <          return
274 <       endif
275 <    
262 >    do i = 1,  EAMList%currentAddition
263 >       number_r   = EAMList%EAMParams(i)%eam_nr
264 >       number_rho = EAMList%EAMParams(i)%eam_nrho
265 >
266 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
267 >            EAMList%EAMParams(i)%eam_rho_r, &
268 >            EAMList%EAMParams(i)%eam_rho_r_pp, &
269 >            0.0E0_DP, 0.0E0_DP, 'N')
270 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
271 >            EAMList%EAMParams(i)%eam_Z_r, &
272 >            EAMList%EAMParams(i)%eam_Z_r_pp, &
273 >            0.0E0_DP, 0.0E0_DP, 'N')
274 >       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
275 >            EAMList%EAMParams(i)%eam_F_rho, &
276 >            EAMList%EAMParams(i)%eam_F_rho_pp, &
277 >            0.0E0_DP, 0.0E0_DP, 'N')
278 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
279 >            EAMList%EAMParams(i)%eam_phi_r, &
280 >            EAMList%EAMParams(i)%eam_phi_r_pp, &
281 >            0.0E0_DP, 0.0E0_DP, 'N')
282 >    enddo
283 >
284 >    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
285 >    !! find the smallest rcut for any eam atype
286 >    !       do i = 2, EAMList%currentAddition
287 >    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
288 >    !       end do
289 >
290 >    !       EAM_rcut = current_rcut_max
291 >    !       EAM_rcut_orig = current_rcut_max
292 >    !       do i = 1, EAMList%currentAddition
293 >    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
294 >    !       end do
295 >    !! Allocate arrays for force calculation
296 >
297 >    call allocateEAM(alloc_stat)
298 >    if (alloc_stat /= 0 ) then
299 >       write(*,*) "allocateEAM failed"
300 >       status = -1
301 >       return
302 >    endif
303 >
304    end subroutine init_EAM_FF
305  
306 < !! routine checks to see if array is allocated, deallocates array if allocated
307 < !! and then creates the array to the required size
306 >  !! routine checks to see if array is allocated, deallocates array if allocated
307 >  !! and then creates the array to the required size
308    subroutine allocateEAM(status)
309      integer, intent(out) :: status
310  
# Line 289 | Line 347 | contains
347         status = -1
348         return
349      end if
350 <    
350 >
351   #ifdef IS_MPI
352  
353      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 326 | Line 384 | contains
384      end if
385  
386  
387 < ! Now do column arrays
387 >    ! Now do column arrays
388  
389      if (allocated(frho_col)) deallocate(frho_col)
390      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
# Line 352 | Line 410 | contains
410         status = -1
411         return
412      end if
413 <  
413 >
414   #endif
415  
416    end subroutine allocateEAM
417  
418 < !! C sets rcut to be the largest cutoff of any atype
419 < !! present in this simulation. Doesn't include all atypes
420 < !! sim knows about, just those in the simulation.
418 >  !! C sets rcut to be the largest cutoff of any atype
419 >  !! present in this simulation. Doesn't include all atypes
420 >  !! sim knows about, just those in the simulation.
421    subroutine setCutoffEAM(rcut, status)
422      real(kind=dp) :: rcut
423      integer :: status
# Line 372 | Line 430 | contains
430  
431  
432    subroutine clean_EAM()
433 <  
434 <   ! clean non-IS_MPI first
433 >
434 >    ! clean non-IS_MPI first
435      frho = 0.0_dp
436      rho  = 0.0_dp
437      dfrhodrho = 0.0_dp
438 < ! clean MPI if needed
438 >    ! clean MPI if needed
439   #ifdef IS_MPI
440      frho_row = 0.0_dp
441      frho_col = 0.0_dp
# Line 401 | Line 459 | contains
459  
460  
461      if (present(stat)) stat = 0
462 <    
462 >
463      allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
464      if (alloc_stat /= 0 ) then
465         if (present(stat)) stat = -1
# Line 452 | Line 510 | contains
510         if (present(stat)) stat = -1
511         return
512      end if
455 –      
513  
514 +
515    end subroutine allocate_EAMType
516  
517  
# Line 461 | Line 519 | contains
519      type (EAMtype), pointer :: thisEAMType
520  
521      ! free Arrays in reverse order of allocation...
522 <    deallocate(thisEAMType%eam_phi_r_pp)      
523 <    deallocate(thisEAMType%eam_rho_r_pp)  
524 <    deallocate(thisEAMType%eam_Z_r_pp)  
525 <    deallocate(thisEAMType%eam_F_rho_pp)  
526 <    deallocate(thisEAMType%eam_phi_r)      
527 <    deallocate(thisEAMType%eam_rho_r)      
528 <    deallocate(thisEAMType%eam_Z_r)  
529 <    deallocate(thisEAMType%eam_F_rho)
530 <    deallocate(thisEAMType%eam_rhovals)
531 <    deallocate(thisEAMType%eam_rvals)
532 <  
522 >    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
523 >    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
524 >    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
525 >    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
526 >    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
527 >    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
528 >    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
529 >    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
530 >    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
531 >    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
532 >
533    end subroutine deallocate_EAMType
534  
535 < !! Calculates rho_r
535 >  !! Calculates rho_r
536    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
537      integer :: atom1,atom2
538      real(kind = dp), dimension(3) :: d
# Line 488 | Line 546 | contains
546      ! we don't use the derivatives, dummy variables
547      real( kind = dp) :: drho,d2rho
548      integer :: eam_err
549 <  
549 >
550      integer :: myid_atom1
551      integer :: myid_atom2
552  
553 < ! check to see if we need to be cleaned at the start of a force loop
496 <      
553 >    ! check to see if we need to be cleaned at the start of a force loop
554  
555  
556  
557 +
558   #ifdef IS_MPI
559      myid_atom1 = atid_Row(atom1)
560      myid_atom2 = atid_Col(atom2)
# Line 516 | Line 574 | contains
574              r, rho_i_at_j,drho,d2rho)
575  
576  
577 <      
577 >
578   #ifdef  IS_MPI
579         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
580   #else
581         rho(atom2) = rho(atom2) + rho_i_at_j
582   #endif
583 < !       write(*,*) atom1,atom2,r,rho_i_at_j
584 <       endif
583 >       !       write(*,*) atom1,atom2,r,rho_i_at_j
584 >    endif
585  
586 <       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
587 <          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
588 <               EAMList%EAMParams(myid_atom2)%eam_rvals, &
589 <               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
590 <               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
591 <               r, rho_j_at_i,drho,d2rho)
592 <
593 <    
594 <      
595 <      
596 < #ifdef  IS_MPI
597 <          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
586 >    if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
587 >       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
588 >            EAMList%EAMParams(myid_atom2)%eam_rvals, &
589 >            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
590 >            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
591 >            r, rho_j_at_i,drho,d2rho)
592 >
593 >
594 >
595 >
596 > #ifdef  IS_MPI
597 >       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
598   #else
599 <          rho(atom1) = rho(atom1) + rho_j_at_i
599 >       rho(atom1) = rho(atom1) + rho_j_at_i
600   #endif
601 <       endif
601 >    endif
602  
603  
604  
# Line 563 | Line 621 | contains
621      integer :: me
622      integer :: n_rho_points
623  
624 <  
624 >
625      cleanme = .true.
626      !! Scatter the electron density from  pre-pair calculation back to local atoms
627   #ifdef IS_MPI
628      call scatter(rho_row,rho,plan_atom_row,eam_err)
629      if (eam_err /= 0 ) then
630 <      write(errMsg,*) " Error scattering rho_row into rho"
631 <      call handleError(RoutineName,errMesg)
632 <   endif      
630 >       write(errMsg,*) " Error scattering rho_row into rho"
631 >       call handleError(RoutineName,errMesg)
632 >    endif
633      call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
634      if (eam_err /= 0 ) then
635 <      write(errMsg,*) " Error scattering rho_col into rho"
636 <      call handleError(RoutineName,errMesg)
637 <   endif
635 >       write(errMsg,*) " Error scattering rho_col into rho"
636 >       call handleError(RoutineName,errMesg)
637 >    endif
638  
639 <      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
639 >    rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
640   #endif
641  
642  
643  
644 < !! Calculate F(rho) and derivative
644 >    !! Calculate F(rho) and derivative
645      do atom = 1, nlocal
646         me = atid(atom)
647         n_rho_points = EAMList%EAMParams(me)%eam_nrho
# Line 613 | Line 671 | contains
671  
672      enddo
673  
616 –  
674  
675 +
676   #ifdef IS_MPI
677      !! communicate f(rho) and derivatives back into row and column arrays
678      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 644 | Line 702 | contains
702      endif
703   #endif
704  
705 <  
705 >
706    end subroutine calc_eam_preforce_Frho
707  
708  
# Line 662 | Line 720 | contains
720      real( kind = dp ), intent(inout), dimension(3) :: fpair
721  
722      logical, intent(in) :: do_pot
723 <    
723 >
724      real( kind = dp ) :: drdx,drdy,drdz
725      real( kind = dp ) :: d2
726      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
# Line 680 | Line 738 | contains
738      integer  :: mytype_atom1
739      integer  :: mytype_atom2
740  
741 < !Local Variables
684 <    
685 <   ! write(*,*) "Frho: ", Frho(atom1)
741 >    !Local Variables
742  
743 +    ! write(*,*) "Frho: ", Frho(atom1)
744 +
745      phab = 0.0E0_DP
746      dvpdr = 0.0E0_DP
747      d2vpdrdr = 0.0E0_DP
# Line 701 | Line 759 | contains
759         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
760         ! get type specific cutoff for atom 2
761         rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
762 <      
762 >
763         drdx = d(1)/rij
764         drdy = d(2)/rij
765         drdz = d(3)/rij
766 <      
766 >
767         if (rij.lt.rci) then
768            call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
769                 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
# Line 727 | Line 785 | contains
785                 EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
786                 EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
787                 rij, rhb,drhb,d2rhb)
788 <          
788 >
789            !! Calculate Phi(r) for atom2.
790            call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
791                 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
# Line 745 | Line 803 | contains
803                 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
804                 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
805         endif
806 <      
806 >
807         if (rij.lt.rcj) then
808            phab = phab + 0.5E0_DP*(rha/rhb)*phb
809            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 755 | Line 813 | contains
813                 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
814                 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
815         endif
816 <      
816 >
817         drhoidr = drha
818         drhojdr = drhb
819  
# Line 770 | Line 828 | contains
828   #else
829         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
830              + dvpdr
831 <      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
831 >       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
832   #endif
833  
834         fx = dudr * drdx
# Line 787 | Line 845 | contains
845         f_Row(1,atom1) = f_Row(1,atom1) + fx
846         f_Row(2,atom1) = f_Row(2,atom1) + fy
847         f_Row(3,atom1) = f_Row(3,atom1) + fz
848 <      
848 >
849         f_Col(1,atom2) = f_Col(1,atom2) - fx
850         f_Col(2,atom2) = f_Col(2,atom2) - fy
851         f_Col(3,atom2) = f_Col(3,atom2) - fz
# Line 800 | Line 858 | contains
858         f(1,atom1) = f(1,atom1) + fx
859         f(2,atom1) = f(2,atom1) + fy
860         f(3,atom1) = f(3,atom1) + fz
861 <      
861 >
862         f(1,atom2) = f(1,atom2) - fx
863         f(2,atom2) = f(2,atom2) - fy
864         f(3,atom2) = f(3,atom2) - fz
865   #endif
866 <      
866 >
867         vpair = vpair + phab
868   #ifdef IS_MPI
869         id1 = AtomRowToGlobal(atom1)
# Line 814 | Line 872 | contains
872         id1 = atom1
873         id2 = atom2
874   #endif
875 <      
875 >
876         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
877 <          
877 >
878            fpair(1) = fpair(1) + fx
879            fpair(2) = fpair(2) + fy
880            fpair(3) = fpair(3) + fz
881 <          
881 >
882         endif
883 <    
883 >
884         if (nmflag) then
885  
886            drhoidr = drha
# Line 836 | Line 894 | contains
894                 d2rhojdrdr*dfrhodrho_row(atom1) + &
895                 drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
896                 drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
897 <              
897 >
898   #else
899  
900            d2 = d2vpdrdr + &
# Line 846 | Line 904 | contains
904                 drhojdr*drhojdr*d2frhodrhodrho(atom1)
905   #endif
906         end if
907 <      
908 <    endif    
907 >
908 >    endif
909    end subroutine do_eam_pair
910  
911  
# Line 862 | Line 920 | contains
920      real( kind = DP ) :: del, h, a, b, c, d
921      integer :: pp_arraySize
922  
923 <
923 >
924      ! this spline code assumes that the x points are equally spaced
925      ! do not attempt to use this code if they are not.
926 <    
927 <    
926 >
927 >
928      ! find the closest point with a value below our own:
929      j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
930  
# Line 878 | Line 936 | contains
936      ! check to make sure we haven't screwed up the calculation of j:
937      if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
938         if (j.ne.nx) then
939 <        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
940 <       call handleError(routineName,errMSG)
939 >          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
940 >          call handleError(routineName,errMSG)
941         endif
942      endif
943  
944      del = xa(j+1) - x
945      h = xa(j+1) - xa(j)
946 <    
946 >
947      a = del / h
948      b = 1.0E0_DP - a
949      c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
950      d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
951 <    
951 >
952      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 –  
953  
954 +    dy = (ya(j+1)-ya(j))/h &
955 +         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
956 +         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
957 +
958 +
959 +    d2y = a*yppa(j) + b*yppa(j+1)
960 +
961 +
962    end subroutine eam_splint
963  
964  
# Line 912 | Line 970 | contains
970      ! if boundary is 'U' the upper derivative is used
971      ! if boundary is 'B' then both derivatives are used
972      ! if boundary is anything else, then both derivatives are assumed to be 0
973 <    
973 >
974      integer :: nx, i, k, max_array_size
975 <    
975 >
976      real( kind = DP ), dimension(:)        :: xa
977      real( kind = DP ), dimension(:)        :: ya
978      real( kind = DP ), dimension(:)        :: yppa
979      real( kind = DP ), dimension(size(xa)) :: u
980      real( kind = DP ) :: yp1,ypn,un,qn,sig,p
981      character(len=*) :: boundary
982 <    
982 >
983      ! make sure the sizes match
984      if ((nx /= size(xa)) .or. (nx /= size(ya))) then
985         call handleWarning("EAM_SPLINE","Array size mismatch")
# Line 936 | Line 994 | contains
994         yppa(1) = 0.0E0_DP
995         u(1)  = 0.0E0_DP
996      endif
997 <    
997 >
998      do i = 2, nx - 1
999         sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1000         p = sig * yppa(i-1) + 2.0E0_DP
# Line 945 | Line 1003 | contains
1003              (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1004              (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1005      enddo
1006 <    
1006 >
1007      if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1008           (boundary.eq.'b').or.(boundary.eq.'B')) then
1009         qn = 0.5E0_DP
# Line 957 | Line 1015 | contains
1015      endif
1016  
1017      yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1018 <    
1018 >
1019      do k = nx-1, 1, -1
1020         yppa(k)=yppa(k)*yppa(k+1)+u(k)
1021      enddo
1022  
1023    end subroutine eam_spline
1024  
967 –
968 –
969 –
1025   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