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

Diff Legend

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