« back to all changes in this revision
Viewing changes to Src/write_subs.F90
- browse files at revision 1070
- compare with another revision
- download diff
- download tarball
- view history from revision 1070
- Committer: Nick Papior
- Date: 2019-02-14 10:00:14 UTC
- Revision ID: nickpapior@gmail.com-20190214100014-pnax8vnrua8363ep
Added 1D and 2D stress tensor calculations
The problem with low-D materials is that stress is "per volume",
but since vacuum is independent in a slab calculation the stress will
be variable depending on the amount of vacuum.
What one *really* should do is to use the volume of the slab along the
vacuum direction so it may be fixed depending on the extend of the slab.
This poses additional problems since how do you decide the slab thickness?
- vdW bond-lengths
- orbital radii
- or 3rd?
This amends output to do the simplest thing but forces users to do post-processing.
The vacuum region is divided out and thus change the units depending on the dimensionality.
- 2D == eV/Ang^2
- 1D == eV/Ang
Then the user can them-selves determine the actual height/area of the slab/chain.
Note that the regular stress tensors are still printed out, but now an additional
block is added.
TODO: these stress's should probably be the ones used for relaxations and var-cell
calculations?
The problem with low-D materials is that stress is "per volume",
but since vacuum is independent in a slab calculation the stress will
be variable depending on the amount of vacuum.
What one *really* should do is to use the volume of the slab along the
vacuum direction so it may be fixed depending on the extend of the slab.
This poses additional problems since how do you decide the slab thickness?
- vdW bond-lengths
- orbital radii
- or 3rd?
This amends output to do the simplest thing but forces users to do post-processing.
The vacuum region is divided out and thus change the units depending on the dimensionality.
- 2D == eV/Ang^2
- 1D == eV/Ang
Then the user can them-selves determine the actual height/area of the slab/chain.
Note that the regular stress tensors are still printed out, but now an additional
block is added.
TODO: these stress's should probably be the ones used for relaxations and var-cell
calculations?
- files modified:
- Src/write_subs.F90
added
removed
Src/write_subs.F90
69
69
70
70
subroutine siesta_write_stress_pressure()
71
71
use m_steps, only: final
72
use siesta_geom, only: nsc
73
74
integer :: nD
75
logical :: periodic(3)
76
77
! TODO, probably we need to use shaper...
78
! However, the shaper returns vectors which are
79
! not parallel to the actual lattice vectors and hence
80
! it becomes difficult to utilize them
81
! What to do for Gamma calculations?
82
periodic(:) = nsc(:) > 1
83
nD = count( periodic )
72
84
73
85
if ( final ) then
74
call write_stress_pressure_final()
86
call write_stress_pressure_final(nD, periodic)
75
87
else
76
call write_stress_pressure_not_final()
88
call write_stress_pressure_not_final(nD, periodic)
77
89
end if
78
90
79
91
end subroutine siesta_write_stress_pressure
80
92
81
subroutine write_stress_pressure_not_final()
93
subroutine write_stress_pressure_not_final(nD, periodic)
82
94
use precision
83
95
use siesta_cml
84
96
use units
89
101
use m_energies, only: FreeE
90
102
use m_stress, only: stress, kin_stress, mstress, tstress
91
103
104
integer, intent(in) :: nD ! dimensionality of the simulation
105
logical, intent(in) :: periodic(3)
92
106
integer :: ix
93
107
94
108
real(dp):: Psol ! Pressure of "solid"
95
109
real(dp):: Press ! Pressure
96
110
real(dp):: ps(3,3) ! Auxiliary array
97
111
98
112
! Stress tensor and pressure:
99
113
100
114
! Write Voigt components of total stress tensor
142
156
' ',stress(:,1)*Ang**3/eV, &
143
157
' ',stress(:,2)*Ang**3/eV, &
144
158
' ',stress(:,3)*Ang**3/eV
159
160
select case ( nD )
161
case ( 1 )
162
call calc_stress_1D(periodic, stress, ps)
163
write(6,'(/,a,3(/,a,3f15.6))') &
164
'siesta: 1D stress tensor (static) (eV/Ang):', &
165
' ',ps(:,1)*Ang/eV, &
166
' ',ps(:,2)*Ang/eV, &
167
' ',ps(:,3)*Ang/eV
168
case ( 2 )
169
call calc_stress_2D(periodic, stress, ps)
170
write(6,'(/,a,3(/,a,3f15.6))') &
171
'siesta: 2D stress tensor (static) (eV/Ang**2):', &
172
' ',ps(:,1)*Ang**2/eV, &
173
' ',ps(:,2)*Ang**2/eV, &
174
' ',ps(:,3)*Ang**2/eV
175
end select
145
176
146
177
Psol = - (stress(1,1) + stress(2,2) + stress(3,3)) / 3.0_dp
147
178
write(6,'(/,a,f20.8,a)') 'siesta: Pressure (static):', Psol/kBar, ' kBar'
157
188
' ',tstress(:,2)*Ang**3/eV, &
158
189
' ',tstress(:,3)*Ang**3/eV
159
190
191
select case ( nD )
192
case ( 1 )
193
call calc_stress_1D(periodic, tstress, ps)
194
write(6,'(/,a,3(/,a,3f15.6))') &
195
'siesta: 1D stress tensor (total) (eV/Ang):', &
196
' ',ps(:,1)*Ang/eV, &
197
' ',ps(:,2)*Ang/eV, &
198
' ',ps(:,3)*Ang/eV
199
case ( 2 )
200
call calc_stress_2D(periodic, tstress, ps)
201
write(6,'(/,a,3(/,a,3f15.6))') &
202
'siesta: 2D stress tensor (total) (eV/Ang**2):', &
203
' ',ps(:,1)*Ang**2/eV, &
204
' ',ps(:,2)*Ang**2/eV, &
205
' ',ps(:,3)*Ang**2/eV
206
end select
207
160
208
Psol = - (tstress(1,1) + tstress(2,2) + tstress(3,3)) / 3.0_dp
161
209
write(6,'(/,a,f20.8,a)') 'siesta: Pressure (total):', Psol/kBar, ' kBar'
162
210
215
263
216
264
end subroutine write_stress_pressure_not_final
217
265
218
subroutine write_stress_pressure_final()
266
subroutine write_stress_pressure_final(nD, periodic)
219
267
use precision
220
268
use siesta_cml
221
269
use units
224
272
use m_iostruct, only: write_struct
225
273
use m_stress, only: stress, mstress, cstress
226
274
275
integer, intent(in) :: nD ! dimensionality of the simulation
276
logical, intent(in) :: periodic(3)
277
227
278
real(dp):: Pmol ! Molecular pressure (discounting Virial term)
228
279
real(dp):: Psol ! Pressure of "solid"
280
real(dp) :: ps(3,3) ! auxiliary array
281
229
282
230
283
! AG: Possible BUG
231
284
! The volume here refers to the "old" cell, and
241
294
dictref='siesta:stress', units='siestaUnits:eV_Ang__3')
242
295
end if
243
296
297
select case ( nD )
298
case ( 1 )
299
call calc_stress_1D(periodic, stress, ps)
300
write(6,'(/,a,3(/,a,3f15.6))') &
301
'siesta: 1D stress tensor (static) (eV/Ang):', &
302
'siesta: ',ps(:,1)*Ang/eV, &
303
'siesta: ',ps(:,2)*Ang/eV, &
304
'siesta: ',ps(:,3)*Ang/eV
305
case ( 2 )
306
call calc_stress_2D(periodic, stress, ps)
307
write(6,'(/,a,3(/,a,3f15.6))') &
308
'siesta: 2D stress tensor (static) (eV/Ang**2):', &
309
'siesta: ',ps(:,1)*Ang**2/eV, &
310
'siesta: ',ps(:,2)*Ang**2/eV, &
311
'siesta: ',ps(:,3)*Ang**2/eV
312
end select
313
244
314
! Print constrained stress tensor if different from unconstrained
245
315
if ( any(cstress /= stress ) ) then
246
316
write(6,'(/,a,3(/,a,3f12.6))') 'siesta: Constrained stress tensor (static) (eV/Ang**3):', &
252
322
dictref='siesta:cstress', &
253
323
units='siestaUnits:eV_Ang__3')
254
324
end if
325
326
select case ( nD )
327
case ( 1 )
328
call calc_stress_1D(periodic, cstress, ps)
329
write(6,'(/,a,3(/,a,3f15.6))') &
330
'siesta: 1D constrained stress tensor (static) (eV/Ang):', &
331
'siesta: ',ps(:,1)*Ang/eV, &
332
'siesta: ',ps(:,2)*Ang/eV, &
333
'siesta: ',ps(:,3)*Ang/eV
334
case ( 2 )
335
call calc_stress_2D(periodic, cstress, ps)
336
write(6,'(/,a,3(/,a,3f15.6))') &
337
'siesta: 2D constrained stress tensor (static) (eV/Ang**2):', &
338
'siesta: ',ps(:,1)*Ang**2/eV, &
339
'siesta: ',ps(:,2)*Ang**2/eV, &
340
'siesta: ',ps(:,3)*Ang**2/eV
341
end select
342
255
343
end if
256
344
257
345
! Find pressure
281
369
end if
282
370
283
371
end subroutine write_stress_pressure_final
372
373
!< Transfer the stress from a chain calculation from eV / Ang ** 3 -> eV / Ang
374
subroutine calc_stress_1D(periodic, stress, stress_1d)
375
use precision
376
! TODO for varcel it isn't necessarily volume_of_some_cell == VOLUME(ucell)
377
use siesta_geom, only: volume_of_some_cell, ucell
378
379
logical, intent(in) :: periodic(3)
380
real(dp), intent(in) :: stress(3,3)
381
real(dp), intent(out) :: stress_1d(3,3)
382
383
real(dp) :: vlen
384
385
stress_1d(:,:) = 0._dp
386
if ( periodic(1) ) then
387
stress_1d(1,1) = stress(1,1)
388
vlen = sqrt( dot_product(ucell(:,1), ucell(:,1)) )
389
else if ( periodic(2) ) then
390
stress_1d(2,2) = stress(2,2)
391
vlen = sqrt( dot_product(ucell(:,2), ucell(:,2)) )
392
else if ( periodic(3) ) then
393
stress_1d(3,3) = stress(3,3)
394
vlen = sqrt( dot_product(ucell(:,3), ucell(:,3)) )
395
else
396
call die('write_subs_pressure: 1D error in periodic determination')
397
end if
398
399
stress_1d(:,:) = stress_1d(:, :) * volume_of_some_cell / vlen
400
401
end subroutine calc_stress_1D
402
403
!< Transfer the stress from a chain calculation from eV / Ang ** 3 -> eV / Ang ** 2
404
subroutine calc_stress_2D(periodic, stress, stress_2d)
405
use precision
406
! TODO for varcel it isn't necessarily volume_of_some_cell == VOLUME(ucell)
407
use siesta_geom, only: volume_of_some_cell, ucell
408
409
logical, intent(in) :: periodic(3)
410
real(dp), intent(in) :: stress(3,3)
411
real(dp), intent(out) :: stress_2d(3,3)
412
real(dp) :: vcross(3)
413
414
stress_2d(:,:) = stress(:,:)
415
if ( periodic(1) .and. periodic(2) ) then
416
stress_2d(:,3) = 0._dp
417
stress_2d(3,:) = 0._dp
418
call cross(ucell(:,1), ucell(:,2), vcross)
419
else if ( periodic(1) .and. periodic(3) ) then
420
stress_2d(:,2) = 0._dp
421
stress_2d(2,:) = 0._dp
422
call cross(ucell(:,1), ucell(:,3), vcross)
423
else if ( periodic(2) .and. periodic(3) ) then
424
stress_2d(:,1) = 0._dp
425
stress_2d(1,:) = 0._dp
426
call cross(ucell(:,2), ucell(:,3), vcross)
427
else
428
call die('write_subs_pressure: 2D error in periodic determination')
429
end if
430
431
stress_2d(:,:) = stress_2d(:,:) * volume_of_some_cell / sqrt( dot_product(vcross, vcross) )
432
433
end subroutine calc_stress_2D
284
434
285
435
end module write_subs_pressure
286
436
Loggerhead is a web-based interface for Breezy
Version: 2.0.1