ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE/libmdtools/calc_eam.F90

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

Diff Legend

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