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

Diff Legend

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