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 2277 by chrisfen, Fri Aug 26 21:30:41 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
# 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  
# Line 110 | Line 154 | contains
154      integer                                :: eam_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", eam_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
144 <    ! EAMAtypes
145 <      
191 >  
192      EAMList%EAMParams(current)%eam_atype    = eam_ident
193      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
194      EAMList%EAMParams(current)%eam_nrho     = eam_nrho
# 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 +    real(kind=dp) :: cutValue
226 +    
227 +    cutValue = eamList%EAMParams(atomID)%eam_rcut
228 +
229 +  end function getEAMCut
230 +
231    subroutine init_EAM_FF(status)
232      integer :: status
233      integer :: i,j
# Line 173 | Line 243 | contains
243         return
244      end if
245  
176 –
177 –       do i = 1, EAMList%currentAddition
246  
247 < ! Build array of r values
247 >    do i = 1, EAMList%currentAddition
248  
249 <          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
249 >       ! Build array of r values
250  
251 <          ! precompute the pair potential and get it into kcal / mol:
252 <          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
253 <          do j = 2, EAMList%EAMParams(i)%eam_nr
254 <             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
251 >       do j = 1,EAMList%EAMParams(i)%eam_nr
252 >          EAMList%EAMParams(i)%eam_rvals(j) = &
253 >               real(j-1,kind=dp)* &
254 >               EAMList%EAMParams(i)%eam_dr
255         end do
256 <      
256 >       ! Build array of rho values
257 >       do j = 1,EAMList%EAMParams(i)%eam_nrho
258 >          EAMList%EAMParams(i)%eam_rhovals(j) = &
259 >               real(j-1,kind=dp)* &
260 >               EAMList%EAMParams(i)%eam_drho
261 >       end do
262 >       ! convert from eV to kcal / mol:
263 >       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
264  
265 <       do i = 1,  EAMList%currentAddition
266 <          number_r   = EAMList%EAMParams(i)%eam_nr
267 <          number_rho = EAMList%EAMParams(i)%eam_nrho
268 <          
269 <          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')
265 >       ! precompute the pair potential and get it into kcal / mol:
266 >       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
267 >       do j = 2, EAMList%EAMParams(i)%eam_nr
268 >          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
269 >          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
270         enddo
271 +    end do
272  
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
273  
274 < !       EAM_rcut = current_rcut_max
275 < !       EAM_rcut_orig = current_rcut_max
276 < !       do i = 1, EAMList%currentAddition
277 < !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
278 < !       end do
279 <       !! Allocate arrays for force calculation
280 <      
281 <       call allocateEAM(alloc_stat)
282 <       if (alloc_stat /= 0 ) then
283 <          write(*,*) "allocateEAM failed"
284 <          status = -1
285 <          return
286 <       endif
287 <    
274 >    do i = 1,  EAMList%currentAddition
275 >       number_r   = EAMList%EAMParams(i)%eam_nr
276 >       number_rho = EAMList%EAMParams(i)%eam_nrho
277 >
278 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
279 >            EAMList%EAMParams(i)%eam_rho_r, &
280 >            EAMList%EAMParams(i)%eam_rho_r_pp, &
281 >            0.0E0_DP, 0.0E0_DP, 'N')
282 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
283 >            EAMList%EAMParams(i)%eam_Z_r, &
284 >            EAMList%EAMParams(i)%eam_Z_r_pp, &
285 >            0.0E0_DP, 0.0E0_DP, 'N')
286 >       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
287 >            EAMList%EAMParams(i)%eam_F_rho, &
288 >            EAMList%EAMParams(i)%eam_F_rho_pp, &
289 >            0.0E0_DP, 0.0E0_DP, 'N')
290 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
291 >            EAMList%EAMParams(i)%eam_phi_r, &
292 >            EAMList%EAMParams(i)%eam_phi_r_pp, &
293 >            0.0E0_DP, 0.0E0_DP, 'N')
294 >    enddo
295 >
296 >    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
297 >    !! find the smallest rcut for any eam atype
298 >    !       do i = 2, EAMList%currentAddition
299 >    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
300 >    !       end do
301 >
302 >    !       EAM_rcut = current_rcut_max
303 >    !       EAM_rcut_orig = current_rcut_max
304 >    !       do i = 1, EAMList%currentAddition
305 >    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
306 >    !       end do
307 >    !! Allocate arrays for force calculation
308 >
309 >    call allocateEAM(alloc_stat)
310 >    if (alloc_stat /= 0 ) then
311 >       write(*,*) "allocateEAM failed"
312 >       status = -1
313 >       return
314 >    endif
315 >
316    end subroutine init_EAM_FF
317  
318 < !! routine checks to see if array is allocated, deallocates array if allocated
319 < !! and then creates the array to the required size
318 >  !! routine checks to see if array is allocated, deallocates array if allocated
319 >  !! and then creates the array to the required size
320    subroutine allocateEAM(status)
321      integer, intent(out) :: status
322  
# Line 289 | Line 359 | contains
359         status = -1
360         return
361      end if
362 <    
362 >
363   #ifdef IS_MPI
364  
365      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 326 | Line 396 | contains
396      end if
397  
398  
399 < ! Now do column arrays
399 >    ! Now do column arrays
400  
401      if (allocated(frho_col)) deallocate(frho_col)
402      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
# Line 352 | Line 422 | contains
422         status = -1
423         return
424      end if
425 <  
425 >
426   #endif
427  
428    end subroutine allocateEAM
429  
430 < !! C sets rcut to be the largest cutoff of any atype
431 < !! present in this simulation. Doesn't include all atypes
432 < !! sim knows about, just those in the simulation.
430 >  !! C sets rcut to be the largest cutoff of any atype
431 >  !! present in this simulation. Doesn't include all atypes
432 >  !! sim knows about, just those in the simulation.
433    subroutine setCutoffEAM(rcut, status)
434      real(kind=dp) :: rcut
435      integer :: status
# Line 372 | Line 442 | contains
442  
443  
444    subroutine clean_EAM()
445 <  
446 <   ! clean non-IS_MPI first
445 >
446 >    ! clean non-IS_MPI first
447      frho = 0.0_dp
448      rho  = 0.0_dp
449      dfrhodrho = 0.0_dp
450 < ! clean MPI if needed
450 >    ! clean MPI if needed
451   #ifdef IS_MPI
452      frho_row = 0.0_dp
453      frho_col = 0.0_dp
# Line 401 | Line 471 | contains
471  
472  
473      if (present(stat)) stat = 0
474 <    
474 >
475      allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
476      if (alloc_stat /= 0 ) then
477         if (present(stat)) stat = -1
# Line 452 | Line 522 | contains
522         if (present(stat)) stat = -1
523         return
524      end if
455 –      
525  
526 +
527    end subroutine allocate_EAMType
528  
529  
# Line 461 | Line 531 | contains
531      type (EAMtype), pointer :: thisEAMType
532  
533      ! free Arrays in reverse order of allocation...
534 <    deallocate(thisEAMType%eam_phi_r_pp)      
535 <    deallocate(thisEAMType%eam_rho_r_pp)  
536 <    deallocate(thisEAMType%eam_Z_r_pp)  
537 <    deallocate(thisEAMType%eam_F_rho_pp)  
538 <    deallocate(thisEAMType%eam_phi_r)      
539 <    deallocate(thisEAMType%eam_rho_r)      
540 <    deallocate(thisEAMType%eam_Z_r)  
541 <    deallocate(thisEAMType%eam_F_rho)
542 <    deallocate(thisEAMType%eam_rhovals)
543 <    deallocate(thisEAMType%eam_rvals)
544 <  
534 >    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
535 >    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
536 >    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
537 >    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
538 >    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
539 >    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
540 >    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
541 >    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
542 >    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
543 >    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
544 >
545    end subroutine deallocate_EAMType
546  
547 < !! Calculates rho_r
547 >  !! Calculates rho_r
548    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
549      integer :: atom1,atom2
550      real(kind = dp), dimension(3) :: d
# Line 488 | Line 558 | contains
558      ! we don't use the derivatives, dummy variables
559      real( kind = dp) :: drho,d2rho
560      integer :: eam_err
561 <  
561 >
562      integer :: myid_atom1
563      integer :: myid_atom2
564  
565 < ! check to see if we need to be cleaned at the start of a force loop
496 <      
565 >    ! check to see if we need to be cleaned at the start of a force loop
566  
567  
568  
569 +
570   #ifdef IS_MPI
571      myid_atom1 = atid_Row(atom1)
572      myid_atom2 = atid_Col(atom2)
# Line 516 | Line 586 | contains
586              r, rho_i_at_j,drho,d2rho)
587  
588  
589 <      
589 >
590   #ifdef  IS_MPI
591         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
592   #else
593         rho(atom2) = rho(atom2) + rho_i_at_j
594   #endif
595 < !       write(*,*) atom1,atom2,r,rho_i_at_j
596 <       endif
527 <
528 <       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
529 <          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
530 <               EAMList%EAMParams(myid_atom2)%eam_rvals, &
531 <               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
532 <               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
533 <               r, rho_j_at_i,drho,d2rho)
595 >       !       write(*,*) atom1,atom2,r,rho_i_at_j
596 >    endif
597  
598 <    
599 <      
600 <      
598 >    if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
599 >       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
600 >            EAMList%EAMParams(myid_atom2)%eam_rvals, &
601 >            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
602 >            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
603 >            r, rho_j_at_i,drho,d2rho)
604 >
605 >
606 >
607 >
608   #ifdef  IS_MPI
609 <          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
609 >       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
610   #else
611 <          rho(atom1) = rho(atom1) + rho_j_at_i
611 >       rho(atom1) = rho(atom1) + rho_j_at_i
612   #endif
613 <       endif
613 >    endif
614  
615  
616  
# Line 563 | Line 633 | contains
633      integer :: me
634      integer :: n_rho_points
635  
636 <  
636 >
637      cleanme = .true.
638      !! Scatter the electron density from  pre-pair calculation back to local atoms
639   #ifdef IS_MPI
640      call scatter(rho_row,rho,plan_atom_row,eam_err)
641      if (eam_err /= 0 ) then
642 <      write(errMsg,*) " Error scattering rho_row into rho"
643 <      call handleError(RoutineName,errMesg)
644 <   endif      
642 >       write(errMsg,*) " Error scattering rho_row into rho"
643 >       call handleError(RoutineName,errMesg)
644 >    endif
645      call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
646      if (eam_err /= 0 ) then
647 <      write(errMsg,*) " Error scattering rho_col into rho"
648 <      call handleError(RoutineName,errMesg)
649 <   endif
647 >       write(errMsg,*) " Error scattering rho_col into rho"
648 >       call handleError(RoutineName,errMesg)
649 >    endif
650  
651 <      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
651 >    rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
652   #endif
653  
654  
655  
656 < !! Calculate F(rho) and derivative
656 >    !! Calculate F(rho) and derivative
657      do atom = 1, nlocal
658         me = atid(atom)
659         n_rho_points = EAMList%EAMParams(me)%eam_nrho
# Line 613 | Line 683 | contains
683  
684      enddo
685  
616 –  
686  
687 +
688   #ifdef IS_MPI
689      !! communicate f(rho) and derivatives back into row and column arrays
690      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 644 | Line 714 | contains
714      endif
715   #endif
716  
717 <  
717 >
718    end subroutine calc_eam_preforce_Frho
719  
720  
# Line 662 | Line 732 | contains
732      real( kind = dp ), intent(inout), dimension(3) :: fpair
733  
734      logical, intent(in) :: do_pot
735 <    
735 >
736      real( kind = dp ) :: drdx,drdy,drdz
737      real( kind = dp ) :: d2
738      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
# Line 680 | Line 750 | contains
750      integer  :: mytype_atom1
751      integer  :: mytype_atom2
752  
753 < !Local Variables
754 <    
755 <   ! write(*,*) "Frho: ", Frho(atom1)
753 >    !Local Variables
754 >
755 >    ! write(*,*) "Frho: ", Frho(atom1)
756  
757      phab = 0.0E0_DP
758      dvpdr = 0.0E0_DP
# Line 701 | Line 771 | contains
771         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
772         ! get type specific cutoff for atom 2
773         rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
774 <      
774 >
775         drdx = d(1)/rij
776         drdy = d(2)/rij
777         drdz = d(3)/rij
778 <      
778 >
779         if (rij.lt.rci) then
780            call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
781                 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
# Line 727 | Line 797 | contains
797                 EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
798                 EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
799                 rij, rhb,drhb,d2rhb)
800 <          
800 >
801            !! Calculate Phi(r) for atom2.
802            call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
803                 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
# Line 745 | Line 815 | contains
815                 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
816                 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
817         endif
818 <      
818 >
819         if (rij.lt.rcj) then
820            phab = phab + 0.5E0_DP*(rha/rhb)*phb
821            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 755 | Line 825 | contains
825                 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
826                 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
827         endif
828 <      
828 >
829         drhoidr = drha
830         drhojdr = drhb
831  
# Line 770 | Line 840 | contains
840   #else
841         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
842              + dvpdr
843 <      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
843 >       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
844   #endif
845  
846         fx = dudr * drdx
# Line 787 | Line 857 | contains
857         f_Row(1,atom1) = f_Row(1,atom1) + fx
858         f_Row(2,atom1) = f_Row(2,atom1) + fy
859         f_Row(3,atom1) = f_Row(3,atom1) + fz
860 <      
860 >
861         f_Col(1,atom2) = f_Col(1,atom2) - fx
862         f_Col(2,atom2) = f_Col(2,atom2) - fy
863         f_Col(3,atom2) = f_Col(3,atom2) - fz
# Line 800 | Line 870 | contains
870         f(1,atom1) = f(1,atom1) + fx
871         f(2,atom1) = f(2,atom1) + fy
872         f(3,atom1) = f(3,atom1) + fz
873 <      
873 >
874         f(1,atom2) = f(1,atom2) - fx
875         f(2,atom2) = f(2,atom2) - fy
876         f(3,atom2) = f(3,atom2) - fz
877   #endif
878 <      
878 >
879         vpair = vpair + phab
880   #ifdef IS_MPI
881         id1 = AtomRowToGlobal(atom1)
# Line 814 | Line 884 | contains
884         id1 = atom1
885         id2 = atom2
886   #endif
887 <      
887 >
888         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
889 <          
889 >
890            fpair(1) = fpair(1) + fx
891            fpair(2) = fpair(2) + fy
892            fpair(3) = fpair(3) + fz
893 <          
893 >
894         endif
895 <    
895 >
896         if (nmflag) then
897  
898            drhoidr = drha
# Line 836 | Line 906 | contains
906                 d2rhojdrdr*dfrhodrho_row(atom1) + &
907                 drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
908                 drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
909 <              
909 >
910   #else
911  
912            d2 = d2vpdrdr + &
# Line 846 | Line 916 | contains
916                 drhojdr*drhojdr*d2frhodrhodrho(atom1)
917   #endif
918         end if
919 <      
920 <    endif    
919 >
920 >    endif
921    end subroutine do_eam_pair
922  
923  
# Line 862 | Line 932 | contains
932      real( kind = DP ) :: del, h, a, b, c, d
933      integer :: pp_arraySize
934  
935 <
935 >
936      ! this spline code assumes that the x points are equally spaced
937      ! do not attempt to use this code if they are not.
938 <    
939 <    
938 >
939 >
940      ! find the closest point with a value below our own:
941      j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
942  
# Line 878 | Line 948 | contains
948      ! check to make sure we haven't screwed up the calculation of j:
949      if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
950         if (j.ne.nx) then
951 <        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
952 <       call handleError(routineName,errMSG)
951 >          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
952 >          call handleError(routineName,errMSG)
953         endif
954      endif
955  
956      del = xa(j+1) - x
957      h = xa(j+1) - xa(j)
958 <    
958 >
959      a = del / h
960      b = 1.0E0_DP - a
961      c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
962      d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
963 <    
963 >
964      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 –  
965  
966 +    dy = (ya(j+1)-ya(j))/h &
967 +         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
968 +         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
969 +
970 +
971 +    d2y = a*yppa(j) + b*yppa(j+1)
972 +
973 +
974    end subroutine eam_splint
975  
976  
# Line 912 | Line 982 | contains
982      ! if boundary is 'U' the upper derivative is used
983      ! if boundary is 'B' then both derivatives are used
984      ! if boundary is anything else, then both derivatives are assumed to be 0
985 <    
985 >
986      integer :: nx, i, k, max_array_size
987 <    
987 >
988      real( kind = DP ), dimension(:)        :: xa
989      real( kind = DP ), dimension(:)        :: ya
990      real( kind = DP ), dimension(:)        :: yppa
991      real( kind = DP ), dimension(size(xa)) :: u
992      real( kind = DP ) :: yp1,ypn,un,qn,sig,p
993      character(len=*) :: boundary
994 <    
994 >
995      ! make sure the sizes match
996      if ((nx /= size(xa)) .or. (nx /= size(ya))) then
997         call handleWarning("EAM_SPLINE","Array size mismatch")
# Line 936 | Line 1006 | contains
1006         yppa(1) = 0.0E0_DP
1007         u(1)  = 0.0E0_DP
1008      endif
1009 <    
1009 >
1010      do i = 2, nx - 1
1011         sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1012         p = sig * yppa(i-1) + 2.0E0_DP
# Line 945 | Line 1015 | contains
1015              (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1016              (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1017      enddo
1018 <    
1018 >
1019      if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1020           (boundary.eq.'b').or.(boundary.eq.'B')) then
1021         qn = 0.5E0_DP
# Line 957 | Line 1027 | contains
1027      endif
1028  
1029      yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1030 <    
1030 >
1031      do k = nx-1, 1, -1
1032         yppa(k)=yppa(k)*yppa(k+1)+u(k)
1033      enddo
1034  
1035    end subroutine eam_spline
1036  
967 –
968 –
969 –
1037   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