-
Notifications
You must be signed in to change notification settings - Fork 0
/
output.f90
441 lines (371 loc) · 12.7 KB
/
output.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
module output
!
! Purpose :
! Contains output routines
!
! Date Author History of Revison
! ==== ====== ==================
! 18.02.2014 Svenja M. Janke Original
! Sascha Kandratsenka
! Dan J. Auerbach
use open_file
use atom_class
use md_init
implicit none
save
integer :: save_counter = 1
logical :: overwrite = .true.
contains
subroutine full_conf(slab, teil, itraj,Eref) !mxt_conf.bin
!
! Purpose:
! Prints out information necessary to continue simulation
! Results in binary file.
!
type(atoms) :: slab, teil
real(8) :: Eref
integer :: ios, itraj
character(len=8) str
character(len=90) filename
write(str,'(I8.8)') save_counter
filename = 'conf/mxt_conf'//str//'.bin'
open (753,file=filename,form="unformatted", status='replace', &
action='write', iostat=ios)
! Here comes the output
write(753) itraj ! Number of trajectory
write(753) step ! time step
write(753) Epot, Eref
write(753) Tsurf ! Surface temperature
! number of species
if (teil%n_atoms .ne. 0) then
write(753) 2
else
write(753) 1
end if
! potential and neighbouring
write(753) pes_name!, pes_nigh
! name, number of atoms, no of fixed atoms
! masses, number of parameters, file name parameters, paramter values
! propagator
write(753) name_l, slab%n_atoms, slab%nofix,&
mass_l, npars_l,key_l
write(753 )pars_l, md_algo_l
write(753) a_lat ! lattice constant
write(753) cell_mat ! Cell matrix
write(753) cell_imat ! inverse cell matrix
write(753) slab%r, slab%v, slab%a, slab%dens
if (teil%n_atoms .ne. 0) then
write(753) name_p, teil%n_atoms, teil%nofix, &
mass_p, npars_p, key_p
write(753) pars_p, md_algo_p
write(753) teil%r, teil%v, teil%a, teil%dens
end if
close(753)
save_counter = save_counter+1
end subroutine full_conf
subroutine out_short(slab, teil,Epot, Eref, itraj, q, rmin_p, col_int, imp, rbounce) !mxt_fin.dat
!
! Purpose :
! Prints out final information at end of trajectory
! Results in mxt_fin-files necessary for traj analysis
!
type(atoms) :: slab, teil
real(8) :: Epot, Ekin_l, Ekin_p, Eref
integer :: itraj, q, i
real(8), dimension(:,:), allocatable :: rmin_p
real(8), dimension(:,:,:), allocatable :: rbounce
integer, dimension(:), allocatable :: col_int, imp
character(len=8) str
character(len=90) filename
write(str,'(I8.8)') itraj
filename = 'traj/mxt_fin'//str//'.dat'
if (overwrite) then
call open_for_write(753,filename)
else
call open_for_append(753,filename)
write(753,*)
write(753,'(A14,f15.5)') 'Time (fs) = ', q*step
end if
Ekin_l = E_kin(slab,mass_l)
Ekin_p = E_kin(teil,mass_p)
write(753,'(A14,f15.5)') 'E_kin_p (eV) = ', Ekin_p
write(753,'(A14,f15.5)') 'E_kin_l (eV) = ', Ekin_l
write(753,'(A14,f15.5)') 'E_pot (eV) = ', Epot
if (md_algo_p == 5) then
write(753,'(A14,f15.5)') 'E_pef (eV) = ', pEfric
elseif (md_algo_p == 3 .or. md_algo_p == 4) then
write(753,'(A14,f15.5)') 'E_phonon(eV) = ', pEfric
end if
write(753,'(A14,f15.5)') 'E_ref (eV) = ', Eref
write(753,'(A14,f15.5)') 'E_total (eV) = ', Epot + Ekin_l + Ekin_p
write(753,'(A8)')'r_p (A):'
write(753,'(3f15.5)') teil%r
write(753,'(A11)')'v_p (A/fs):'
write(753,'(3e15.5)') teil%v
if (.not. overwrite) then
write(753,'(A12)')'r_min_p (A):'
write(753,'(3f15.5)') rmin_p
write(753,'(A21)')'Time at surface (fs):'
write(753,'(1000f15.5)') col_int*step
write(753,'(A21)')'Number of bounces:'
write(753,'(1000I10)') imp
write(753,'(A21)')'Bounce sites (A):'
do i=1,teil%n_atoms
if (imp(i)< 5) then
write(753,'(15f15.5)') rbounce(1:imp(i),:,i)
else
write(753,'(15f15.5)') rbounce(1:5,:,i)
end if
end do
end if
close(753)
overwrite = .false.
end subroutine out_short
subroutine out_detail(output_info, n, itraj,Eref) !mxt_trj.dat
!
! Purpose :
! Prints out a lot of trajectory information along trajectory
! mxt_trj-files to follow what's happening during trajectory
!
real(8), dimension(:,:), allocatable :: output_info
integer :: n, i
integer :: itraj
character(len=8) str
character(len=90) filename
real(8) :: Eref
write(str,'(I8.8)') itraj
filename = 'traj/mxt_trj'//str//'.dat'
call open_for_write(753,filename)
write(753,'(1000A15)') 'time(fs)', 'E_pot(eV)', 'E_kin_l(eV)','E_kin_p(eV)',&
'E_total(eV)', 'dens(A^-3)', 'r_p(A)', 'v_p(A/fs)', &
'eta_p (1/fs)'
do i = 1, n
write(753,'(10000e15.5)') i*wstep(2)*step, output_info(:,i)
end do
write(753,'(A15,e15.5)') 'E_ref (eV) = ', Eref
close(753)
end subroutine out_detail
subroutine out_all(slab, teil, itraj) !mxt_conf.xyz
!
! Purpose:
! Prints out all the geometries along the trajectory.
! saved as .xyz-file
! Be careful. This option will take up a lot of disc-space.
!
type(atoms) :: slab, teil
integer :: ios, itraj,i, n
character(len=8) str
character(len=90) filename
real(8), allocatable, dimension(:,:) :: xdatx
write(str,'(I8.8)') itraj
filename = 'conf/mxt_conf'//str//'.xyz'
n=slab%n_atoms+teil%n_atoms
allocate(xdatx(3,n+1))
xdatx=0.0d0
xdatx(1,1:slab%n_atoms)= slab%r(1,:)
xdatx(2,1:slab%n_atoms)= slab%r(2,:)
xdatx(3,1:slab%n_atoms)= slab%r(3,:)
xdatx(1,slab%n_atoms+1:n)=teil%r(1,:)
xdatx(2,slab%n_atoms+1:n)=teil%r(2,:)
xdatx(3,slab%n_atoms+1:n)=teil%r(3,:)
! Get all atoms into unit-cell
xdatx = matmul(cell_imat,xdatx)
do i = 1, n
! if (slab%r(1,i) > 1) slab%r(1,i) = slab%r(1,i) - int(slab%r(1,i))
! if (slab%r(2,i) > 1) slab%r(2,i) = slab%r(2,i) - int(slab%r(2,i))
! if (slab%r(1,i) < 0.0d0) slab%r(1,i) = slab%r(1,i) - int(slab%r(1,i))+1
! if (slab%r(2,i) < 0.0d0) slab%r(2,i) = slab%r(2,i) - int(slab%r(2,i))+1
if (xdatx(1,i) > 1.100d0) xdatx(1,i) = xdatx(1,i) - Anint(xdatx(1,i))
if (xdatx(2,i) > 1.1d0) xdatx(2,i) = xdatx(2,i) - Anint(xdatx(2,i))
if (xdatx(1,i) < -1.1d0) xdatx(1,i) = xdatx(1,i) - Anint(xdatx(1,i))+1
if (xdatx(2,i) < -1.1d0) xdatx(2,i) = xdatx(2,i) - Anint(xdatx(2,i))+1
end do
xdatx = matmul(cell_mat,xdatx)
! open (753,file=filename, status='replace', &
! action='write', iostat=ios)
if (overwrite) then
open (753,file=filename, status='replace', &
action='write', iostat=ios)
else
call open_for_append(753,filename)
end if
! Here comes the output
write(753,*) slab%n_atoms+teil%n_atoms
write(753,*)
do i = 1, slab%n_atoms
write(753,*) 'Au ', xdatx(:,i)
end do
do i = slab%n_atoms +1 , slab%n_atoms+teil%n_atoms
write(753,*) 'H ', xdatx(:,i)
end do
save_counter = save_counter+1
close(753)
overwrite = .false.
deallocate(xdatx)
end subroutine out_all
subroutine out_poscar(slab,teil,Epot, Eref, itraj) !mxt_anneal.POSCAR
!
! Purpose :
! Prints out poscar-file for final state
!
type(atoms) :: slab, teil
real(8) :: Epot, Ekin_l, Ekin_p, Eref
integer :: itraj
character(len=8) str
character(len=90) filename
write(str,'(I8.8)') itraj
filename = 'traj/mxt_anneal'//str//'.POSCAR'
call open_for_write(753,filename)
Ekin_l = E_kin(slab,mass_l)
Ekin_p = E_kin(teil,mass_p)
write(753,'(5(A6,f10.5),A5)') 'Ek_p=', Ekin_p, 'Ek_l=', Ekin_l, &
'Epot=', Epot, 'Eref=', Eref, 'Etot=', Epot + Ekin_l + Ekin_p, '(eV)'
write(753,*) 1.0
write(753,'(3f18.8)') cell_mat
if(teil%n_atoms>0) then
write(753,*) slab%n_atoms,teil%n_atoms
else
write(753,*) slab%n_atoms
end if
write(753,*) 'Cartesian'
write(753,'(3f18.8)') slab%r
if(teil%n_atoms>0) write(753,'(3f18.8)') teil%r
close(753)
end subroutine out_poscar
function sartre(itraj)
!
! Purpose:
! Checks if trajectory has been calculated before
! If following trajectory has been calculated, too, we want to omit
! calculating this one by setting sartre to true.
!
logical :: sartre
integer :: itraj
logical :: exists
character(len=8) str
character(len=80) filename, filenamen
write(str,'(I8.8)') itraj
sartre = .false.
if (wstep(1) == -1) filename='traj/mxt_fin'//str//'.dat'
if (wstep(1) == 0) filename='traj/mxt_trj'//str//'.dat'
inquire(file=filename ,exist=exists)
if (exists) then
write(str,'(I8.8)') itraj+1
if (wstep(1) == -1) filenamen='traj/mxt_fin'//str//'.dat'
if (wstep(1) == 0) filenamen='traj/mxt_trj'//str//'.dat'
inquire(file=filenamen ,exist=exists)
if (exists) then
sartre = .true.
else
sartre = .false.
filename = 'rm -f '//trim(filename)
call system(filename)
end if
else
sartre = .false.
end if
end function sartre
subroutine out_posvel(slab, teil, itraj, Eref) !mxt_rv.dat
!
! Purpose:
! Prints out all the geometries and velocities along the trajectory.
! Be careful. This option will take up a lot of space.
!
type(atoms) :: slab, teil
real(8) :: Eref
integer :: ios, itraj
character(len=8) str
character(len=90) filename
write(str,'(I8.8)') save_counter
filename = 'conf/mxt_rv'//str//'.dat'
open (753,file=filename, status='replace', &
action='write', iostat=ios)
! Here comes the output
write(753,*) itraj ! Number of trajectory
write(753,*) Epot, Eref
write(753,*) Tsurf ! Surface temperature
! ! number of species
! if (teil%n_atoms .ne. 0) then
! write(753,*) 2
! else
! write(753,*) 1
! end if
! name, number of atoms, no of fixed atoms
! masses, number of parameters, file name parameters, paramter values
! propagator
write(753,*) name_l, slab%n_atoms, slab%nofix,&
mass_l, npars_l,key_l
write(753,*) pars_l, md_algo_l
write(753,*) a_lat ! lattice constant
write(753,*) cell_mat ! Cell matrix
write(753,*) slab%r, slab%v, slab%dens !slab%a
if (teil%n_atoms .ne. 0) then
write(753,*) name_p, teil%n_atoms, teil%nofix, &
mass_p, npars_p, key_p
write(753,*) pars_p, md_algo_p
write(753,*) teil%r, teil%v,teil%dens !, teil%a,
end if
close(753)
!filename = 'gzip '//filename
!call system(filename)
save_counter = save_counter+1
end subroutine out_posvel
subroutine out_pdb(slab, teil, q) !mxt_conf.pdb
!
! Purpose:
! Prints out all the geometries along the trajectory.
! saved as .pdb-file
! Be careful. This option will take up a lot of disc-space.
!
type(atoms) :: slab, teil
integer :: ios, q,i, n
character(len=8) str
character(len=90) filename
! real(8), allocatable, dimension(:,:) :: xdatx
write(str,'(I8)') q
filename = 'conf/mxt_conf'//str//'.pdb'
n=slab%n_atoms+teil%n_atoms
if (overwrite) then
open (753,file=filename, status='replace', &
action='write', iostat=ios)
else
call open_for_append(753,filename)
end if
! Here comes the output
write(753,'(A6,3f9.3,3f7.2, A11)') 'CRYST1', cell_mat(1,1),&
cell_mat(1,1),cell_mat(3,3), 90.00,90.00,120.00,'f m 3 m 4'
do i=1,slab%n_atoms
write(753,'(A6,I5, A1, A4,A1,A3,A2,I4,A4,3f8.3,2f6.2,A10,A2,A2)') &
'ATOM ', i, &
' ',&
'Au',&
' ',&
'UNK',&
' ',&
1,&
' ',&
slab%r(1,i),slab%r(2,i),slab%r(3,i),1.00,0.00,&
' ',&
'Au',&
' '
end do
do i=1,teil%n_atoms
write(753,'(A6,I5, A1, A4,A1,A3,A2,I4,A4,3f8.3,2f6.2,A10,A2,A2)') &
'ATOM ', i, &
' ',&
'H',&
' ',&
'UNK',&
' ',&
1,&
' ',&
teil%r(1,i),teil%r(2,i),teil%r(3,i),1.00,0.00,&
' ',&
' H',&
' '
end do
save_counter = save_counter+1
close(753)
overwrite = .false.
end subroutine out_pdb
end module output