ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_eam.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_eam.F90 (file contents):
Revision 497 by chuckv, Mon Apr 14 21:16:37 2003 UTC vs.
Revision 631 by chuckv, Thu Jul 17 19:25:51 2003 UTC

# Line 1 | Line 1
1 < module calc_eam
2 <  use definitions, ONLY : DP
1 > module eam
2 >  use definitions, ONLY : DP,default_error
3 >  use simulation
4    use force_globals
5 +  use status
6 +  use atype_module
7   #ifdef MPI
8    use mpiSimulation
9   #endif
10 +  implicit none
11    PRIVATE
12  
13  
14    logical, save :: EAM_FF_initialized = .false.
15    integer, save :: EAM_Mixing_Policy
16 <  integer, save :: EAM_rcut
16 >  real(kind = dp), save :: EAM_rcut
17  
18 +  character(len = 200) :: errMsg
19 +  character(len=*), parameter :: RoutineName =  "EAM MODULE"
20 +  logical :: cleanme = .true.
21  
22  
23 +
24 +  type, private :: EAMtype
25 +     integer           :: eam_atype      
26 +     real( kind = DP ) :: eam_dr          
27 +     integer           :: eam_nr          
28 +     integer           :: eam_nrho          
29 +     real( kind = DP ) :: eam_lattice        
30 +     real( kind = DP ) :: eam_drho      
31 +     real( kind = DP ) :: eam_rcut    
32 +     integer           :: eam_atype_map
33 +    
34 +     real( kind = DP ), pointer, dimension(:) :: eam_rvals        => null()
35 +     real( kind = DP ), pointer, dimension(:) :: eam_rhovals      => null()
36 +     real( kind = DP ), pointer, dimension(:) :: eam_F_rho        => null()
37 +     real( kind = DP ), pointer, dimension(:) :: eam_Z_r          => null()
38 +     real( kind = DP ), pointer, dimension(:) :: eam_rho_r        => null()
39 +     real( kind = DP ), pointer, dimension(:) :: eam_phi_r        => null()
40 +     real( kind = DP ), pointer, dimension(:) :: eam_F_rho_pp     => null()
41 +     real( kind = DP ), pointer, dimension(:) :: eam_Z_r_pp       => null()
42 +     real( kind = DP ), pointer, dimension(:) :: eam_rho_r_pp     => null()
43 +     real( kind = DP ), pointer, dimension(:) :: eam_phi_r_pp     => null()
44 +  end type EAMtype
45  
17 –  !! standard eam stuff  
18 –  integer                                        :: n_eam_atypes
19 –  integer,           allocatable, dimension(:)   :: eam_atype    
20 –  real( kind = DP ), allocatable, dimension(:)   :: eam_dr      
21 –  integer, allocatable, dimension(:)   :: eam_nr      
22 –  integer, allocatable, dimension(:)   :: eam_nrho      
23 –  real( kind = DP ), allocatable, dimension(:)   :: eam_lattice      
24 –  real( kind = DP ), allocatable, dimension(:)   :: eam_drho    
25 –  integer          , allocatable, dimension(:)   :: eam_atype_map
26 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rvals    
27 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rhovals  
28 –  real( kind = DP ), allocatable, dimension(:)   :: eam_rcut    
29 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_F_rho    
30 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_Z_r      
31 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rho_r    
32 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_phi_r    
33 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_F_rho_pp
34 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_Z_r_pp  
35 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rho_r_pp
36 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_phi_r_pp
46  
47 +  !! Arrays for derivatives used in force calculation
48 +  real( kind = dp), dimension(:), allocatable :: frho
49 +  real( kind = dp), dimension(:), allocatable :: rho
50  
51 <  public :: init_EAM_FF
52 <  public :: EAM_new_rcut
41 <  public :: do_EAM_pair
51 >  real( kind = dp), dimension(:), allocatable :: dfrhodrho
52 > !  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
53  
54  
55 + !! Arrays for MPI storage
56 + #ifdef MPI
57 +  real( kind = dp), dimension(:), allocatable :: dfrhodrho_col
58 +  real( kind = dp), dimension(:), allocatable :: dfrhodrho_row
59 +  real( kind = dp), dimension(:), allocatable :: frho_row
60 +  real( kind = dp), dimension(:), allocatable :: frho_col
61 +  real( kind = dp), dimension(:), allocatable :: rho_row
62 +  real( kind = dp), dimension(:), allocatable :: rho_col
63  
64 < contains
46 <  subroutine init_EAM_FF()
64 > #endif
65  
66 +  type, private :: EAMTypeList
67 +     integer           :: n_eam_types = 0
68 +     integer           :: currentAddition = 0
69 +    
70 +     type (EAMtype), pointer  :: EAMParams(:) => null()
71 +  end type EAMTypeList
72  
73  
74 <    character(len=80) :: eam_pot_file
51 <    integer :: i, j,  max_size, prev_max_size
52 <    integer :: number_rho, number_r
53 <    integer :: eam_unit
54 <    integer :: this_error
55 <    character(len=300) :: msg
56 <    integer, external :: nfiles
57 <    !for mpi
74 >  type (eamTypeList) :: EAMList
75  
76 +  !! standard eam stuff  
77  
60 – #ifdef MPI
61 –    if (node == 0)  &  
62 –         n_eam_atypes = nfiles(trim(eam_pot_dir)//char(0))  
78  
79 <    call mpi_bcast(n_eam_atypes,1,mpi_integer,0,mpi_comm_world,mpi_err)
80 <    if (n_eam_atypes == -1) then
81 <       call error("INITIALIZE_EAM","NO EAM potentials found!")
82 <    endif
68 <    write(msg,'(a5,i4,a12,i5,a14)') 'Node: ',node,' reading ...', &
69 <         n_eam_atypes, ' eam atom types'
70 <    call info('INITIALIZE_EAM', trim(msg))
71 < #else
72 <    n_eam_atypes = nfiles(trim(eam_pot_dir)//char(0))
73 <    if (n_eam_atypes == -1) then
74 <       call error("INITIALIZE_EAM","NO EAM potentials found!")
75 <    endif
79 >  public :: init_EAM_FF
80 > !  public :: EAM_new_rcut
81 > !  public :: do_EAM_pair
82 >  public :: newEAMtype
83  
84 <    write(msg,'(a12,i5,a14)') ' Reading ...', &
78 <         n_eam_atypes, ' eam atom types'
79 <    call info('INITIALIZE_EAM', trim(msg))
80 < #endif
84 >  
85  
86 + contains
87  
83 –    call allocate_eam_atype(n_eam_atypes)
88  
89 +  subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
90 +       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
91 +       eam_ident,status)
92 +    real (kind = dp )                      :: lattice_constant
93 +    integer                                :: eam_nrho
94 +    real (kind = dp )                      :: eam_drho
95 +    integer                                :: eam_nr
96 +    real (kind = dp )                      :: eam_dr
97 +    real (kind = dp )                      :: rcut
98 +    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r
99 +    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
100 +    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
101 +    integer                                :: eam_ident
102 +    integer                                :: status
103  
104 +    integer                                :: nAtypes
105 +    integer                                :: maxVals
106 +    integer                                :: alloc_stat
107 +    integer                                :: current
108 +    integer,pointer                        :: Matchlist(:) => null()
109 +    status = 0
110  
111 <    !! get largest number of data points for any potential
112 < #ifdef MPI
113 <    if (node == 0) then
114 < #endif
115 <       prev_max_size = 0
116 <       do i = 1, n_eam_atypes
117 <          call getfilename(i, eam_pot_file)
118 <          max_size = max(get_eam_sizes( &
119 <               trim(eam_pot_dir) // '/' // eam_pot_file), &
96 <               prev_max_size)
97 <          prev_max_size = max_size
98 <       end do
99 < #ifdef MPI
111 >    !! Assume that atypes has already been set and get the total number of types in atypes
112 >    !! Also assume that every member of atypes is a EAM model.
113 >  
114 >
115 >    ! check to see if this is the first time into
116 >    if (.not.associated(EAMList%EAMParams)) then
117 >       call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
118 >       EAMList%n_eam_types = nAtypes
119 >       allocate(EAMList%EAMParams(nAtypes))
120      end if
121  
122 +    EAMList%currentAddition = EAMList%currentAddition + 1
123 +    current = EAMList%currentAddition
124 +    
125  
126 <    call mpi_bcast(max_size,1,mpi_integer,0,mpi_comm_world,mpi_err)
127 < #endif
126 >    !call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
127 >    if (alloc_stat /= 0) then
128 >       status = -1
129 >       return
130 >    end if
131  
106 –    call allocate_eam_module(n_eam_atypes,max_size)
107 –    allocate(eam_atype_map(get_max_atype()))
132  
133 < #ifdef MPI
134 <    if (node == 0) then
135 < #endif
136 <       do i = 1, n_eam_atypes
137 <          call getfilename(i, eam_pot_file)
138 <          call read_eam_pot(i,trim(eam_pot_dir) // '/' // eam_pot_file, &
139 <               this_error)
133 >    ! this is a possible bug, we assume a correspondence between the vector atypes and
134 >    ! EAMAtypes
135 >      
136 >    EAMList%EAMParams(current)%eam_atype    = eam_ident
137 >    EAMList%EAMParams(current)%eam_lattice  = lattice_constant
138 >    EAMList%EAMParams(current)%eam_nrho     = eam_nrho
139 >    EAMList%EAMParams(current)%eam_drho     = eam_drho
140 >    EAMList%EAMParams(current)%eam_nr       = eam_nr
141 >    EAMList%EAMParams(current)%eam_rcut     = rcut
142 >    EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
143 >    EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
144 >    EAMList%EAMParams(current)%eam_F_rho    = eam_F_rho
145  
146 <          do j = 1, eam_nr(i)
118 <             eam_rvals(j,i) = dble(j-1)*eam_dr(i)
119 <          enddo
146 >  end subroutine newEAMtype
147  
121 –          do j = 1, eam_nrho(i)
122 –             eam_rhovals(j,i) = dble(j-1)*eam_drho(i)
123 –          enddo
148  
149 +
150 +  subroutine init_EAM_FF(status)
151 +    integer :: status
152 +    integer :: i,j
153 +    real(kind=dp) :: current_rcut_max
154 +    integer :: alloc_stat
155 +    integer :: number_r, number_rho
156 +
157 +       do i = 1, EAMList%currentAddition
158 +  
159 +          EAMList%EAMParams(i)%eam_rvals(1:EAMList%EAMParams(i)%eam_nr) = &
160 +               real(EAMList%EAMParams(i)%eam_nr,kind=dp)* &
161 +               EAMList%EAMParams(i)%eam_dr
162 +
163 +          EAMList%EAMParams(i)%eam_rhovals(1:EAMList%EAMParams(i)%eam_nrho) = &
164 +               real(EAMList%EAMParams(i)%eam_nrho,kind=dp)* &
165 +               EAMList%EAMParams(i)%eam_drho
166 +
167            ! convert from eV to kcal / mol:
168 <          do j = 1, eam_nrho(i)
127 <             eam_F_rho(j,i) = eam_F_rho(j,i)*23.06054E0_DP
128 <          enddo
168 >          EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
169  
170            ! precompute the pair potential and get it into kcal / mol:
171 <          eam_phi_r(1,i) = 0.0E0_DP
172 <          do j = 2, eam_nr(i)
173 <             eam_phi_r(j,i) = (eam_Z_r(j,i)**2)/eam_rvals(j,i)
174 <             eam_phi_r(j,i) = eam_phi_r(j,i)*331.999296E0_DP
171 >          EAMList%EAMParams(i)%eam_phi_r = 0.0E0_DP
172 >          do j = 2, EAMList%EAMParams(i)%eam_nr
173 >             EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
174 >             EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
175            enddo
176 +       end do
177  
178 +       do i = 1,  EAMList%currentAddition
179 +          number_r   = EAMList%EAMParams(i)%eam_nr
180 +          number_rho = EAMList%EAMParams(i)%eam_nrho
181 +          
182 +          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
183 +               EAMList%EAMParams(i)%eam_rho_r, &
184 +               EAMList%EAMParams(i)%eam_rho_r_pp, &
185 +               0.0E0_DP, 0.0E0_DP, 'N')
186 +          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
187 +               EAMList%EAMParams(i)%eam_Z_r, &
188 +               EAMList%EAMParams(i)%eam_Z_r_pp, &
189 +               0.0E0_DP, 0.0E0_DP, 'N')
190 +          call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
191 +               EAMList%EAMParams(i)%eam_F_rho, &
192 +               EAMList%EAMParams(i)%eam_F_rho_pp, &
193 +               0.0E0_DP, 0.0E0_DP, 'N')
194 +          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
195 +               EAMList%EAMParams(i)%eam_phi_r, &
196 +               EAMList%EAMParams(i)%eam_phi_r_pp, &
197 +               0.0E0_DP, 0.0E0_DP, 'N')
198 +       enddo
199 +      
200 +       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
201 +       !! find the smallest rcut
202 +       do i = 2, EAMList%currentAddition
203 +          current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
204         end do
138 – #ifdef MPI
139 –       call info('INITIALIZE_EAM','NODE 0: Distributing spline arrays')
140 –    endif
205  
206 <    call mpi_bcast(this_error,n_eam_atypes,mpi_integer,0, &
207 <         mpi_comm_world,mpi_err)
208 <    if (this_error /= 0) then
209 <       call error('INITIALIZE_EAM',"Cannot read eam files")
146 <    endif
206 >       EAM_rcut = current_rcut_max
207 > !       do i = 1, EAMList%currentAddition
208 > !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
209 > !       end do
210  
148 –    call mpi_bcast(eam_atype,n_eam_atypes,mpi_integer,0, &
149 –         mpi_comm_world,mpi_err)
211  
212 <    !! distribute values to cluster......
213 <    call mpi_bcast(eam_nr,n_eam_atypes,mpi_integer,&
214 <         0,mpi_comm_world,mpi_err)
215 <    call mpi_bcast(eam_nrho,n_eam_atypes,mpi_integer,&
216 <         0,mpi_comm_world,mpi_err)
217 <    call mpi_bcast(eam_rvals,n_eam_atypes*max_size,mpi_double_precision, &
157 <         0,mpi_comm_world,mpi_err)
158 <    call mpi_bcast(eam_rcut,n_eam_atypes,mpi_double_precision, &
159 <         0,mpi_comm_world,mpi_err)
160 <    call mpi_bcast(eam_rhovals,n_eam_atypes*max_size,mpi_double_precision, &
161 <         0,mpi_comm_world,mpi_err)
212 >       !! Allocate arrays for force calculation
213 >          call allocateEAM(alloc_stat)
214 >          if (alloc_stat /= 0 ) then
215 >             status = -1
216 >             return
217 >          endif
218  
219 <    !! distribute arrays
164 <    call mpi_bcast(eam_rho_r,n_eam_atypes*max_size,mpi_double_precision, &
165 <         0,mpi_comm_world,mpi_err)
166 <    call mpi_bcast(eam_Z_r,n_eam_atypes*max_size,mpi_double_precision, &
167 <         0,mpi_comm_world,mpi_err)
168 <    call mpi_bcast(eam_F_rho,n_eam_atypes*max_size,mpi_double_precision, &
169 <         0,mpi_comm_world,mpi_err)
170 <    call mpi_bcast(eam_phi_r,n_eam_atypes*max_size,mpi_double_precision, &
171 <         0,mpi_comm_world,mpi_err)
219 >  end subroutine init_EAM_FF
220  
221 < #endif
222 <    call info('INITIALIZE_EAM', 'creating splines')
223 <
224 <    do i = 1, n_eam_atypes
177 <       number_r = eam_nr(i)
178 <       number_rho = eam_nrho(i)
179 <
180 <       call eam_spline(i, number_r, eam_rvals, eam_rho_r, eam_rho_r_pp, &
181 <            0.0E0_DP, 0.0E0_DP, 'N')
182 <       call eam_spline(i, number_r, eam_rvals, eam_Z_r, eam_Z_r_pp, &
183 <            0.0E0_DP, 0.0E0_DP, 'N')
184 <       call eam_spline(i, number_rho, eam_rhovals, eam_F_rho, eam_F_rho_pp, &
185 <            0.0E0_DP, 0.0E0_DP, 'N')
186 <       call eam_spline(i, number_r, eam_rvals, eam_phi_r, eam_phi_r_pp, &
187 <            0.0E0_DP, 0.0E0_DP, 'N')
188 <    enddo
189 <
190 <    do i = 1, n_eam_atypes
191 <       eam_atype_map(eam_atype(i)) = i
192 <    end do
221 > !! routine checks to see if array is allocated, deallocates array if allocated
222 > !! and then creates the array to the required size
223 >  subroutine allocateEAM(status)
224 >    integer, intent(out) :: status
225  
226 +    integer :: nlocal
227 + #ifdef MPI
228 +    integer :: nrow
229 +    integer :: ncol
230 + #endif
231 +    integer :: alloc_stat
232  
233  
234 <    call info('INITIALIZE_EAM','Done creating splines')
234 >    nlocal = getNlocal()
235  
236 <    return
237 <  end subroutine initialize_eam
236 > #ifdef MPI
237 >    nrow = getNrow(plan_row)
238 >    ncol = getNcol(plan_col)
239 > #endif
240  
241 +    if (allocated(frho)) deallocate(frho)
242 +    allocate(frho(nlocal),stat=alloc_stat)
243 +    if (alloc_stat /= 0) then
244 +       status = -1
245 +       return
246 +    end if
247 +    if (allocated(rho)) deallocate(rho)
248 +    allocate(rho(nlocal),stat=alloc_stat)
249 +    if (alloc_stat /= 0) then
250 +       status = -1
251 +       return
252 +    end if
253 +    if (allocated(dfrhodrho)) deallocate(dfrhodrho)
254 +    allocate(dfrhodrho(nlocal),stat=alloc_stat)
255 +    if (alloc_stat /= 0) then
256 +       status = -1
257 +       return
258 +    end if
259 +    
260 + #ifdef MPI
261  
262 <  subroutine allocate_eam_atype(n_size_atype)
263 <    integer, intent(in) :: n_size_atype
262 >    if (allocated(frho_row)) deallocate(frho_row)
263 >    allocate(frho_row(nrow),stat=alloc_stat)
264 >    if (alloc_stat /= 0) then
265 >       status = -1
266 >       return
267 >    end if
268 >    if (allocated(rho_row)) deallocate(rho_row)
269 >    allocate(rho_row(nrow),stat=alloc_stat)
270 >    if (alloc_stat /= 0) then
271 >       status = -1
272 >       return
273 >    end if
274 >    if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
275 >    allocate(dfrhodrho_row(nrow),stat=alloc_stat)
276 >    if (alloc_stat /= 0) then
277 >       status = -1
278 >       return
279 >    end if
280  
281 <    allocate(eam_atype(n_size_atype))  
206 <    allocate(eam_drho(n_size_atype))    
207 <    allocate(eam_dr(n_size_atype))      
208 <    allocate(eam_nr(n_size_atype))
209 <    allocate(eam_nrho(n_size_atype))
210 <    allocate(eam_lattice(n_size_atype))      
211 <    allocate(eam_rcut(n_size_atype))  
281 > ! Now do column arrays
282  
283 <  end subroutine allocate_eam_atype
283 >    if (allocated(frho_col)) deallocate(frho_col)
284 >    allocate(frho_col(ncol),stat=alloc_stat)
285 >    if (alloc_stat /= 0) then
286 >       status = -1
287 >       return
288 >    end if
289 >    if (allocated(rho_col)) deallocate(rho_col)
290 >    allocate(rho_col(ncol),stat=alloc_stat)
291 >    if (alloc_stat /= 0) then
292 >       status = -1
293 >       return
294 >    end if
295 >    if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
296 >    allocate(dfrhodrho_col(ncol),stat=alloc_stat)
297 >    if (alloc_stat /= 0) then
298 >       status = -1
299 >       return
300 >    end if
301  
302 <  subroutine allocate_eam_module(n_size_atype,n_eam_points)
216 <    integer, intent(in) :: n_eam_points
217 <    integer, intent(in) :: n_size_atype
302 > #endif
303  
304 <    allocate(eam_rvals(n_eam_points,n_size_atype))  
220 <    allocate(eam_rhovals(n_eam_points,n_size_atype))  
221 <    allocate(eam_F_rho(n_eam_points,n_size_atype))  
222 <    allocate(eam_Z_r(n_eam_points,n_size_atype))        
223 <    allocate(eam_rho_r(n_eam_points,n_size_atype))      
224 <    allocate(eam_phi_r(n_eam_points,n_size_atype))      
225 <    allocate(eam_F_rho_pp(n_eam_points,n_size_atype))  
226 <    allocate(eam_Z_r_pp(n_eam_points,n_size_atype))  
227 <    allocate(eam_rho_r_pp(n_eam_points,n_size_atype))  
228 <    allocate(eam_phi_r_pp(n_eam_points,n_size_atype))  
304 >  end subroutine allocateEAM
305  
230 –  end subroutine allocate_eam_module
306  
307 <  subroutine deallocate_eam_module()
307 >  subroutine clean_EAM()
308  
309 <    deallocate(eam_atype)  
310 <    deallocate(eam_drho)    
311 <    deallocate(eam_dr)      
312 <    deallocate(eam_nr)
313 <    deallocate(eam_nrho)
314 <    deallocate(eam_lattice)      
315 <    deallocate(eam_atype_map)
316 <    deallocate(eam_rvals)  
317 <    deallocate(eam_rhovals)  
318 <    deallocate(eam_rcut)  
319 <    deallocate(eam_Z_r)        
320 <    deallocate(eam_rho_r)      
321 <    deallocate(eam_phi_r)      
322 <    deallocate(eam_F_rho_pp)  
248 <    deallocate(eam_Z_r_pp)  
249 <    deallocate(eam_rho_r_pp)  
250 <    deallocate(eam_phi_r_pp)  
309 > ! clean non-MPI first
310 >    frho = 0.0_dp
311 >    rho  = 0.0_dp
312 >    dfrhodrho = 0.0_dp
313 > ! clean MPI if needed
314 > #ifdef MPI
315 >    frho_row = 0.0_dp
316 >    frho_col = 0.0_dp
317 >    rho_row  = 0.0_dp
318 >    rho_col  = 0.0_dp
319 >    dfrhodrho_row = 0.0_dp
320 >    dfrhodrho_col = 0.0_dp
321 > #endif
322 >  end subroutine clean_EAM
323  
252 –  end subroutine deallocate_eam_module
324  
325  
326 +  subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
327 +    integer, intent(in)          :: eam_n_rho
328 +    integer, intent(in)          :: eam_n_r
329 +    type (EAMType), pointer      :: thisEAMType
330 +    integer, optional   :: stat
331 +    integer             :: alloc_stat
332 +
333 +
334 +
335 +    if (present(stat)) stat = 0
336 +    
337 +    allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
338 +    if (alloc_stat /= 0 ) then
339 +       if (present(stat)) stat = -1
340 +       return
341 +    end if
342 +    allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat)  
343 +    if (alloc_stat /= 0 ) then
344 +       if (present(stat)) stat = -1
345 +       return
346 +    end if
347 +    allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)  
348 +    if (alloc_stat /= 0 ) then
349 +       if (present(stat)) stat = -1
350 +       return
351 +    end if
352 +    allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)        
353 +    if (alloc_stat /= 0 ) then
354 +       if (present(stat)) stat = -1
355 +       return
356 +    end if
357 +    allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat)      
358 +    if (alloc_stat /= 0 ) then
359 +       if (present(stat)) stat = -1
360 +       return
361 +    end if
362 +    allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat)      
363 +    if (alloc_stat /= 0 ) then
364 +       if (present(stat)) stat = -1
365 +       return
366 +    end if
367 +    allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat)  
368 +    if (alloc_stat /= 0 ) then
369 +       if (present(stat)) stat = -1
370 +       return
371 +    end if
372 +    allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)  
373 +    if (alloc_stat /= 0 ) then
374 +       if (present(stat)) stat = -1
375 +       return
376 +    end if
377 +    allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)  
378 +    if (alloc_stat /= 0 ) then
379 +       if (present(stat)) stat = -1
380 +       return
381 +    end if
382 +    allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
383 +    if (alloc_stat /= 0 ) then
384 +       if (present(stat)) stat = -1
385 +       return
386 +    end if
387 +      
388 +
389 +  end subroutine allocate_EAMType
390 +
391 +
392 +  subroutine deallocate_EAMType(thisEAMType)
393 +    type (EAMtype), pointer :: thisEAMType
394 +
395 +    ! free Arrays in reverse order of allocation...
396 +    deallocate(thisEAMType%eam_phi_r_pp)      
397 +    deallocate(thisEAMType%eam_rho_r_pp)  
398 +    deallocate(thisEAMType%eam_Z_r_pp)  
399 +    deallocate(thisEAMType%eam_F_rho_pp)  
400 +    deallocate(thisEAMType%eam_phi_r)      
401 +    deallocate(thisEAMType%eam_rho_r)      
402 +    deallocate(thisEAMType%eam_Z_r)  
403 +    deallocate(thisEAMType%eam_F_rho)
404 +    deallocate(thisEAMType%eam_rhovals)
405 +    deallocate(thisEAMType%eam_rvals)
406 +  
407 +  end subroutine deallocate_EAMType
408 +
409 + !! Calculates rho_r
410 +  subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
411 +    integer :: atom1,atom2
412 +    real(kind = dp), dimension(3) :: d
413 +    real(kind = dp), intent(inout)               :: r
414 +    real(kind = dp), intent(inout)               :: rijsq
415 +    ! value of electron density rho do to atom i at atom j
416 +    real(kind = dp) :: rho_i_at_j
417 +    ! value of electron density rho do to atom j at atom i
418 +    real(kind = dp) :: rho_j_at_i
419 +
420 +    ! we don't use the derivatives, dummy variables
421 +    real( kind = dp) :: drho,d2rho
422 +    integer :: eam_err
423 +  
424 +    if (cleanme) call clean_EAM
425 +    cleanme = .false.
426 +
427 +    call  calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1)
428 +
429 + #ifdef  MPI
430 +    rho_col(atom2) = rho_col(atom2) + rho_i_at_j
431 + #else
432 +    rho(atom2) = rho(atom2) + rho_i_at_j
433 + #endif
434 +
435 +    call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2)
436 +
437 + #ifdef  MPI
438 +    rho_row(atom1) = rho_row(atom1) + rho_j_at_i
439 + #else
440 +    rho(atom1) = rho(atom1) + rho_j_at_i
441 + #endif
442 +  end subroutine calc_eam_prepair_rho
443 +
444 +  !! Calculate the functional F(rho) for all atoms
445 +  subroutine calc_eam_prepair_Frho(nlocal,pot)
446 +    integer :: nlocal
447 +    real(kind=dp) :: pot
448 +    integer :: i,j
449 +    real(kind=dp) :: U,U1,U2
450 +    integer :: atype1
451 +    ! reset clean forces to be true at top of calc rho.
452 +    cleanme = .true.
453 +
454 + !! Scatter the electron density in pre-pair
455 + #ifdef MPI
456 +    call scatter(rho_row,rho,plan_row,eam_err)
457 +    if (eam_err /= 0 ) then
458 +      write(errMsg,*) " Error scattering rho_row into rho"
459 +      call handleError(RoutineName,errMesg)
460 +   endif      
461 +    call scatter(rho_col,rho,plan_col,eam_err)
462 +    if (eam_err /= 0 ) then
463 +      write(errMsg,*) " Error scattering rho_col into rho"
464 +      call handleError(RoutineName,errMesg)
465 +   endif
466 + #endif
467 +
468 +   do i = 1, nlocal
469 +      call calc_eam_frho(rho(i),u,u1,u2,atype1)
470 +      frho(i) = u
471 +      dfrhodrho(i) = u1
472 + !      d2frhodrhodrho(i) = u2
473 +      pot = pot + u
474 +   enddo
475 +
476 + #ifdef MPI
477 +    !! communicate f(rho) and derivatives
478 +    call gather(frho,frho_row,plan_row, eam_err)
479 +    if (eam_err /=  0) then
480 +       call handleError("cal_eam_forces()","MPI gather frho_row failure")
481 +    endif
482 +    call gather(dfrhodrho,dfrhodrho_row,plan_row, eam_err)
483 +    if (eam_err /=  0) then
484 +       call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure")
485 +    endif
486 +    call gather(frho,frho_col,plan_col, eam_err)
487 +    if (eam_err /=  0) then
488 +       call handleError("cal_eam_forces()","MPI gather frho_col failure")
489 +    endif
490 +    call gather(dfrhodrho,dfrhodrho_col,plan_col, eam_err)
491 +    if (eam_err /=  0) then
492 +       call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
493 +    endif
494 +
495 +    if (nmflag) then
496 +       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row)
497 +       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col)
498 +    endif
499 + #endif
500 +
501 +
502 +  end subroutine calc_eam_prepair_Frho
503 +
504 +
505 +
506 +
507    subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
508   !Arguments    
509      integer, intent(in) ::  atom1, atom2
# Line 261 | Line 513 | contains
513      real( kind = dp ), intent(in), dimension(3) :: d
514      logical, intent(in) :: do_pot, do_stress
515  
516 +    real( kind = dp ) :: drdx,drdy,drdz
517 +    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
518 +    real( kind = dp ) :: rha,drha,d2rha, dpha
519 +    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
520 +    real( kind = dp ) :: dudr
521 +    real( kind = dp ) :: rci,rcj
522 +    real( kind = dp ) :: drhoidr,drhojdr
523 +    real( kind = dp ) :: d2rhoidrdr
524 +    real( kind = dp ) :: d2rhojdrdr
525 +    real( kind = dp ) :: Fx,Fy,Fz
526 +    real( kind = dp ) :: r,d2pha,phb,d2phb
527 +    integer :: id1,id2
528 +    integer  :: atype1,atype2
529 +
530 +
531   !Local Variables
532      
533  
534 +    
535 +    phab = 0.0E0_DP
536 +    dvpdr = 0.0E0_DP
537 +    d2vpdrdr = 0.0E0_DP
538 +      
539 +
540      if (rij .lt. EAM_rcut) then
541 + #ifdef IS_MPI
542 + !!!!! FIX ME
543 +       atype1 = atid_row(atom1)
544 + #else
545 +       atype1 = atid(atom1)
546 + #endif
547  
548 <       r = dsqrt(rijsq)
549 <       efr(1,j) = -rxij
550 <       efr(2,j) = -ryij
551 <       efr(3,j) = -rzij
552 <
274 <
548 >       drdx = d(1)/rij
549 >       drdy = d(2)/rij
550 >       drdz = d(3)/rij
551 >      
552 >      
553         call calc_eam_rho(r, rha, drha, d2rha, atype1)
554         call calc_eam_phi(r, pha, dpha, d2pha, atype1)
555 <       rci = eam_rcut(eam_atype_map(atype1))
556 < #ifdef MPI
557 <       atype2 = ident_col(j)
555 >  !     rci = eam_rcut(eam_atype_map(atom1))
556 > #ifdef IS_MPI
557 >       atype2 = atid_col(atom2)
558   #else
559 <       atype2 = ident(j)
559 >       atype2 = atid(atom2)
560   #endif
561 <
561 >      
562         call calc_eam_rho(r, rhb, drhb, d2rhb, atype2)
563         call calc_eam_phi(r, phb, dphb, d2phb, atype2)
564 <       rcj = eam_rcut(eam_atype_map(atype2))
564 >  !     rcj = eam_rcut(eam_atype_map(atype2))
565  
288 –       phab = 0.0E0_DP
289 –       dvpdr = 0.0E0_DP
290 –       d2vpdrdr = 0.0E0_DP
291 –
566         if (r.lt.rci) then
567            phab = phab + 0.5E0_DP*(rhb/rha)*pha
568            dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
# Line 298 | Line 572 | contains
572                 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
573                 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
574         endif
575 +      
576  
302 –
577         if (r.lt.rcj) then
578            phab = phab + 0.5E0_DP*(rha/rhb)*phb
579            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 309 | Line 583 | contains
583                 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
584                 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
585         endif
586 <
313 <
314 < #ifdef MPI
315 <
316 <       e_row(i) = e_row(i) + phab*0.5
317 <       e_col(i) = e_col(i) + phab*0.5
318 < #else
319 <       if (do_pot) pot = pot + phab
320 < #endif
321 <
586 >      
587         drhoidr = drha
588         drhojdr = drhb
589  
590         d2rhoidrdr = d2rha
591         d2rhojdrdr = d2rhb
592 +
593 +
594   #ifdef MPI
595 <       dudr = drhojdr*dfrhodrho_row(i)+drhoidr*dfrhodrho_col(j) &
595 >       dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
596              + dvpdr
597 <
331 <       if (nmflag) then
332 <          d2 = d2vpdrdr + &
333 <               d2rhoidrdr*dfrhodrho_col(j) + &
334 <               d2rhojdrdr*dfrhodrho_row(i) + &
335 <               drhoidr*drhoidr*d2frhodrhodrho_col(j) + &
336 <               drhojdr*drhojdr*d2frhodrhodrho_row(i)
337 <       endif
597 >    
598   #else
599 <       dudr = drhojdr*dfrhodrho(i)+drhoidr*dfrhodrho(j) &
599 >       dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
600              + dvpdr
601  
342 –       d2 = d2vpdrdr + &
343 –            d2rhoidrdr*dfrhodrho(j) + &
344 –            d2rhojdrdr*dfrhodrho(i) + &
345 –            drhoidr*drhoidr*d2frhodrhodrho(j) + &
346 –            drhojdr*drhojdr*d2frhodrhodrho(i)
602   #endif
603  
604 +       fx = dudr * drdx
605 +       fy = dudr * drdy
606 +       fz = dudr * drdz
607  
350 –       do dim = 1, 3                        
608  
352 –          drdx1 = efr(dim,j) / r
353 –          ftmp = dudr * drdx1
354 –
609   #ifdef MPI
610 <          f_col(dim,j) = f_col(dim,j) - ftmp
611 <          f_row(dim,i) = f_row(dim,i) + ftmp
612 < #else
613 <          f(dim,j) = f(dim,j) - ftmp
360 <          f(dim,i) = f(dim,i) + ftmp
361 < #endif
610 >       if (do_pot) then
611 >          pot_Row(atom1) = pot_Row(atom1) + phab*0.5
612 >          pot_Col(atom2) = pot_Col(atom2) + phab*0.5
613 >       end if
614  
615 <          if (nmflag) then
616 <             idim = 3 * (i-1) + dim
617 <             jdim = 3 * (j-1) + dim
615 >       f_Row(1,atom1) = f_Row(1,atom1) + fx
616 >       f_Row(2,atom1) = f_Row(2,atom1) + fy
617 >       f_Row(3,atom1) = f_Row(3,atom1) + fz
618  
619 <             do dim2 = 1, 3
619 >       f_Col(1,atom2) = f_Col(1,atom2) - fx
620 >       f_Col(2,atom2) = f_Col(2,atom2) - fy
621 >       f_Col(3,atom2) = f_Col(3,atom2) - fz
622 > #else
623 >       if(do_pot) pot = pot + phab
624  
625 <                kt1 = d2 * efr(dim,j) * efr(dim2,j)/r/r
626 <                kt2 = - dudr * efr(dim,j) * efr(dim2,j)/r/r/r
625 >       f(1,atom1) = f(1,atom1) + fx
626 >       f(2,atom1) = f(2,atom1) + fy
627 >       f(3,atom1) = f(3,atom1) + fz
628  
629 <                if (dim.eq.dim2) then
630 <                   kt3 = dudr / r
631 <                else
632 <                   kt3 = 0.0E0_DP
376 <                endif
629 >       f(1,atom2) = f(1,atom2) - fx
630 >       f(2,atom2) = f(2,atom2) - fy
631 >       f(3,atom2) = f(3,atom2) - fz
632 > #endif
633  
378 –                ! The factor of 2 below is to compensate for
379 –                ! overcounting.
380 –                ! Mass weighting is done separately...
634  
382 –                ktmp = (kt1+kt2+kt3)/2.0E0_DP
383 –                idim2 = 3 * (i-1) + dim2
384 –                jdim2 = 3 * (j-1) + dim2
635  
386 –                d(idim,  idim2) = d(idim,idim2)  + ktmp
387 –                d(idim2, idim) = d(idim2,idim)   + ktmp
636  
637 <                d(idim,  jdim2) = d(idim,jdim2)  - ktmp
390 <                d(idim2, jdim) = d(idim2,jdim)   - ktmp
637 >       if (do_stress) then
638  
639 <                d(jdim,  idim2) = d(jdim,idim2)  - ktmp
640 <                d(jdim2, idim) = d(jdim2,idim)   - ktmp
639 > #ifdef MPI
640 >          id1 = tagRow(atom1)
641 >          id2 = tagColumn(atom2)
642 > #else
643 >          id1 = atom1
644 >          id2 = atom2
645 > #endif
646  
647 <                d(jdim,  jdim2) = d(jdim,jdim2)  + ktmp
648 <                d(jdim2, jdim) = d(jdim2,jdim)   + ktmp
647 >          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
648 >            
649 >             tau_Temp(1) = tau_Temp(1) + fx * d(1)
650 >             tau_Temp(2) = tau_Temp(2) + fx * d(2)
651 >             tau_Temp(3) = tau_Temp(3) + fx * d(3)
652 >             tau_Temp(4) = tau_Temp(4) + fy * d(1)
653 >             tau_Temp(5) = tau_Temp(5) + fy * d(2)
654 >             tau_Temp(6) = tau_Temp(6) + fy * d(3)
655 >             tau_Temp(7) = tau_Temp(7) + fz * d(1)
656 >             tau_Temp(8) = tau_Temp(8) + fz * d(2)
657 >             tau_Temp(9) = tau_Temp(9) + fz * d(3)
658 >             virial_Temp = virial_Temp + &
659 >                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
660  
398 –             enddo
661            endif
662 <       enddo
663 <
662 >       endif
663 >      
664      endif
403 – enddo
404 – endif
665  
406 – enddo
666  
667 + end subroutine calc_eam_pair
668  
669 + !!$subroutine calc_eam_rho(r, rho, drho, d2rho, atype)
670 + !!$
671 + !!$  !  include 'headers/sizes.h'
672 + !!$
673 + !!$
674 + !!$integer atype, etype, number_r
675 + !!$real( kind = DP )  :: r, rho, drho, d2rho
676 + !!$integer :: i
677 + !!$
678 + !!$
679 + !!$etype = eam_atype_map(atype)
680 + !!$
681 + !!$if (r.lt.eam_rcut(etype)) then
682 + !!$number_r = eam_nr(etype)
683 + !!$call eam_splint(etype, number_r, eam_rvals, eam_rho_r, &
684 + !!$   eam_rho_r_pp, r, rho, drho, d2rho)
685 + !!$else
686 + !!$rho = 0.0E0_DP
687 + !!$drho = 0.0E0_DP
688 + !!$d2rho = 0.0E0_DP
689 + !!$endif
690 + !!$
691 + !!$return
692 + !!$end subroutine calc_eam_rho
693 + !!$
694 + !!$subroutine calc_eam_frho(dens, u, u1, u2, atype)
695 + !!$
696 + !!$  ! include 'headers/sizes.h'
697 + !!$
698 + !!$integer atype, etype, number_rho
699 + !!$real( kind = DP ) :: dens, u, u1, u2
700 + !!$real( kind = DP ) :: rho_vals
701 + !!$
702 + !!$etype = eam_atype_map(atype)
703 + !!$number_rho = eam_nrho(etype)
704 + !!$if (dens.lt.eam_rhovals(number_rho, etype)) then
705 + !!$call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
706 + !!$   eam_f_rho_pp, dens, u, u1, u2)
707 + !!$else
708 + !!$rho_vals = eam_rhovals(number_rho,etype)
709 + !!$call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
710 + !!$   eam_f_rho_pp, rho_vals, u, u1, u2)
711 + !!$endif
712 + !!$
713 + !!$return
714 + !!$end subroutine calc_eam_frho
715 + !!$
716 + !!$subroutine calc_eam_phi(r, phi, dphi, d2phi, atype)
717 + !!$
718 + !!$
719 + !!$
720 + !!$
721 + !!$integer atype, etype, number_r
722 + !!$real( kind = DP ) :: r, phi, dphi, d2phi
723 + !!$
724 + !!$etype = eam_atype_map(atype)
725 + !!$
726 + !!$if (r.lt.eam_rcut(etype)) then
727 + !!$number_r = eam_nr(etype)
728 + !!$call eam_splint(etype, number_r, eam_rvals, eam_phi_r, &
729 + !!$   eam_phi_r_pp, r, phi, dphi, d2phi)
730 + !!$else
731 + !!$phi = 0.0E0_DP
732 + !!$dphi = 0.0E0_DP
733 + !!$d2phi = 0.0E0_DP
734 + !!$endif
735 + !!$
736 + !!$return
737 + !!$end subroutine calc_eam_phi
738  
739  
740 < end subroutine calc_eam_pair
740 >  subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
741  
413 – subroutine calc_eam_rho(r, rho, drho, d2rho, atype)
742  
743 <  !  include 'headers/sizes.h'
743 >    integer :: atype, nx, j
744 >    real( kind = DP ), dimension(:) :: xa
745 >    real( kind = DP ), dimension(:) :: ya
746 >    real( kind = DP ), dimension(:) :: yppa
747 >    real( kind = DP ) :: x, y, dy, d2y
748 >    real( kind = DP ) :: del, h, a, b, c, d
749  
750  
418 – integer atype, etype, number_r
419 – real( kind = DP )  :: r, rho, drho, d2rho
420 – integer :: i
751  
752  
753 < etype = eam_atype_map(atype)
753 >    
754 >    ! this spline code assumes that the x points are equally spaced
755 >    ! do not attempt to use this code if they are not.
756 >    
757 >    
758 >    ! find the closest point with a value below our own:
759 >    j = FLOOR(dble(nx-1) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
760  
761 < if (r.lt.eam_rcut(etype)) then
762 < number_r = eam_nr(etype)
763 < call eam_splint(etype, number_r, eam_rvals, eam_rho_r, &
764 <   eam_rho_r_pp, r, rho, drho, d2rho)
765 < else
766 < rho = 0.0E0_DP
767 < drho = 0.0E0_DP
768 < d2rho = 0.0E0_DP
769 < endif
761 >    ! check to make sure we're inside the spline range:
762 >    if ((j.gt.nx).or.(j.lt.1)) then
763 >       write(default_error,*) "EAM_splint: x is outside bounds of spline"
764 >    endif
765 >    ! check to make sure we haven't screwed up the calculation of j:
766 >    if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
767 >       if (j.ne.nx) then
768 >        write(default_error,*) "EAM_splint: x is outside bounding range"
769 >       endif
770 >    endif
771  
772 < return
773 < end subroutine calc_eam_rho
772 >    del = xa(j+1) - x
773 >    h = xa(j+1) - xa(j)
774 >    
775 >    a = del / h
776 >    b = 1.0E0_DP - a
777 >    c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
778 >    d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
779 >    
780 >    y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
781 >    
782 >    dy = (ya(j+1)-ya(j))/h &
783 >         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
784 >         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
785 >    
786 >    d2y = a*yppa(j) + b*yppa(j+1)
787  
788 < subroutine calc_eam_frho(dens, u, u1, u2, atype)
788 >  end subroutine eam_splint
789  
790 <  ! include 'headers/sizes.h'
790 >  subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
791  
792 < integer atype, etype, number_rho
443 < real( kind = DP ) :: dens, u, u1, u2
444 < real( kind = DP ) :: rho_vals
792 >  
793  
446 – etype = eam_atype_map(atype)
447 – number_rho = eam_nrho(etype)
448 – if (dens.lt.eam_rhovals(number_rho, etype)) then
449 – call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
450 –   eam_f_rho_pp, dens, u, u1, u2)
451 – else
452 – rho_vals = eam_rhovals(number_rho,etype)
453 – call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
454 –   eam_f_rho_pp, rho_vals, u, u1, u2)
455 – endif
794  
795 < return
796 < end subroutine calc_eam_frho
795 >    ! yp1 and ypn are the first derivatives of y at the two endpoints
796 >    ! if boundary is 'L' the lower derivative is used
797 >    ! if boundary is 'U' the upper derivative is used
798 >    ! if boundary is 'B' then both derivatives are used
799 >    ! if boundary is anything else, then both derivatives are assumed to be 0
800 >    
801 >    integer :: nx, i, k, max_array_size
802 >    
803 >    real( kind = DP ), dimension(:)        :: xa
804 >    real( kind = DP ), dimension(:)        :: ya
805 >    real( kind = DP ), dimension(:)        :: yppa
806 >    real( kind = DP ), dimension(size(xa)) :: u
807 >    real( kind = DP ) :: yp1,ypn,un,qn,sig,p
808 >    character boundary
809 >    
810 >    
811 >    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
812 >         (boundary.eq.'b').or.(boundary.eq.'B')) then
813 >       yppa(1) = -0.5E0_DP
814 >       u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
815 >            ya(1))/(xa(2)-xa(1))-yp1)
816 >    else
817 >       yppa(1) = 0.0E0_DP
818 >       u(1)  = 0.0E0_DP
819 >    endif
820 >    
821 >    do i = 2, nx - 1
822 >       sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
823 >       p = sig * yppa(i-1) + 2.0E0_DP
824 >       yppa(i) = (sig - 1.0E0_DP) / p
825 >       u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
826 >            (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
827 >            (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
828 >    enddo
829 >    
830 >    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
831 >         (boundary.eq.'b').or.(boundary.eq.'B')) then
832 >       qn = 0.5E0_DP
833 >       un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
834 >            (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
835 >    else
836 >       qn = 0.0E0_DP
837 >       un = 0.0E0_DP
838 >    endif
839  
840 < subroutine calc_eam_phi(r, phi, dphi, d2phi, atype)
840 >    yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
841 >    
842 >    do k = nx-1, 1, -1
843 >       yppa(k)=yppa(k)*yppa(k+1)+u(k)
844 >    enddo
845  
846 +  end subroutine eam_spline
847  
848  
849  
465 – integer atype, etype, number_r
466 – real( kind = DP ) :: r, phi, dphi, d2phi
850  
851 < etype = eam_atype_map(atype)
469 <
470 < if (r.lt.eam_rcut(etype)) then
471 < number_r = eam_nr(etype)
472 < call eam_splint(etype, number_r, eam_rvals, eam_phi_r, &
473 <   eam_phi_r_pp, r, phi, dphi, d2phi)
474 < else
475 < phi = 0.0E0_DP
476 < dphi = 0.0E0_DP
477 < d2phi = 0.0E0_DP
478 < endif
479 <
480 < return
481 < end subroutine calc_eam_phi
482 <
483 <
484 < subroutine eam_splint(atype, nx, xa, ya, yppa, x, y, dy, d2y)
485 <
486 <  !  include 'headers/sizes.h'
487 <
488 < real( kind = DP ), dimension(:,:) :: xa
489 < real( kind = DP ), dimension(:,:) :: ya
490 < real( kind = DP ), dimension(:,:) :: yppa
491 < real( kind = DP ) :: x, y, dy, d2y
492 < real( kind = DP ) :: del, h, a, b, c, d
493 <
494 <
495 < integer atype, nx, j
496 <
497 <
498 < ! this spline code assumes that the x points are equally spaced
499 < ! do not attempt to use this code if they are not.
500 <
501 <
502 < ! find the closest point with a value below our own:
503 < j = FLOOR(dble(nx-1) * (x - xa(1,atype)) / (xa(nx,atype) - xa(1,atype))) + 1
504 <
505 < ! check to make sure we're inside the spline range:
506 < if ((j.gt.nx).or.(j.lt.1)) call error('eam_splint', &
507 < 'x is outside bounds of spline')
508 <
509 < ! check to make sure we haven't screwed up the calculation of j:
510 < if ((x.lt.xa(j,atype)).or.(x.gt.xa(j+1,atype))) then
511 < if (j.ne.nx) then
512 < call error('eam_splint', &
513 <      'x is outside bounding range')
514 < endif
515 < endif
516 <
517 < del = xa(j+1,atype) - x
518 < h = xa(j+1,atype) - xa(j,atype)
519 <
520 < a = del / h
521 < b = 1.0E0_DP - a
522 < c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
523 < d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
524 <
525 < y = a*ya(j,atype) + b*ya(j+1,atype) + c*yppa(j,atype) + d*yppa(j+1,atype)
526 <
527 < dy = (ya(j+1,atype)-ya(j,atype))/h &
528 < - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j,atype)/6.0E0_DP &
529 < + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1,atype)/6.0E0_DP
530 <
531 < d2y = a*yppa(j,atype) + b*yppa(j+1,atype)
532 <
533 < return
534 < end subroutine eam_splint
535 <
536 < subroutine eam_spline(atype, nx, xa, ya, yppa, yp1, ypn, boundary)
537 <
538 <  !  include 'headers/sizes.h'
539 <
540 <
541 <  ! yp1 and ypn are the first derivatives of y at the two endpoints
542 <  ! if boundary is 'L' the lower derivative is used
543 <  ! if boundary is 'U' the upper derivative is used
544 <  ! if boundary is 'B' then both derivatives are used
545 <  ! if boundary is anything else, then both derivatives are assumed to be 0
546 <
547 < integer nx, i, k, atype, max_array_size
548 <
549 < real( kind = DP ), dimension(:,:) :: xa
550 < real( kind = DP ), dimension(:,:) :: ya
551 < real( kind = DP ), dimension(:,:) :: yppa
552 < real( kind = DP ), allocatable, dimension(:) :: u
553 < real( kind = DP ) :: yp1,ypn,un,qn,sig,p
554 < character boundary
555 <
556 < max_array_size = size(xa,1)
557 < allocate(u(max_array_size))
558 <
559 <
560 < if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
561 < (boundary.eq.'b').or.(boundary.eq.'B')) then
562 < yppa(1, atype) = -0.5E0_DP
563 < u(1) = (3.0E0_DP/(xa(2,atype)-xa(1,atype)))*((ya(2,atype)-&
564 <   ya(1,atype))/(xa(2,atype)-xa(1,atype))-yp1)
565 < else
566 < yppa(1,atype) = 0.0E0_DP
567 < u(1)  = 0.0E0_DP
568 < endif
569 <
570 < do i = 2, nx - 1
571 < sig = (xa(i,atype) - xa(i-1,atype)) / (xa(i+1,atype) - xa(i-1,atype))
572 < p = sig * yppa(i-1,atype) + 2.0E0_DP
573 < yppa(i,atype) = (sig - 1.0E0_DP) / p
574 < u(i) = (6.0E0_DP*((ya(i+1,atype)-ya(i,atype))/(xa(i+1,atype)-xa(i,atype)) - &
575 <   (ya(i,atype)-ya(i-1,atype))/(xa(i,atype)-xa(i-1,atype)))/ &
576 <   (xa(i+1,atype)-xa(i-1,atype)) - sig * u(i-1))/p
577 < enddo
578 <
579 < if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
580 < (boundary.eq.'b').or.(boundary.eq.'B')) then
581 < qn = 0.5E0_DP
582 < un = (3.0E0_DP/(xa(nx,atype)-xa(nx-1,atype)))* &
583 <   (ypn-(ya(nx,atype)-ya(nx-1,atype))/(xa(nx,atype)-xa(nx-1,atype)))
584 < else
585 < qn = 0.0E0_DP
586 < un = 0.0E0_DP
587 < endif
588 <
589 < yppa(nx,atype)=(un-qn*u(nx-1))/(qn*yppa(nx-1,atype)+1.0E0_DP)
590 <
591 < do k = nx-1, 1, -1
592 < yppa(k,atype)=yppa(k,atype)*yppa(k+1,atype)+u(k)
593 < enddo
594 <
595 < deallocate(u)
596 < return
597 < end subroutine eam_spline
598 <
599 <
600 <
601 <
602 < end module calc_eam
851 > end module eam

Diff Legend

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