Branch data Line data Source code
1 : : !!****m* ABINIT/m_ddb_internalstr
2 : : !! NAME
3 : : !! m_ddb_internalstr
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (XW)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_ddb_internalstr
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : use m_crystal
28 : : use m_ddb
29 : :
30 : : use m_fstrings, only : itoa, sjoin
31 : : use m_dynmat, only : asria_corr
32 : :
33 : : implicit none
34 : :
35 : : private
36 : : !!***
37 : :
38 : : public :: ddb_internalstr
39 : : !!***
40 : :
41 : : contains
42 : : !!***
43 : :
44 : : !!****f* ABINIT/ddb_internalstr
45 : : !!
46 : : !! NAME
47 : : !! ddb_internalstr
48 : : !!
49 : : !! FUNCTION
50 : : !! Get the insternal strain tensors,both force response and displacement response ones.
51 : : !!
52 : : !! INPUTS
53 : : !! asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file.
54 : : !! blkval(2,3,mpert,3,mpert,nblok)=
55 : : !! second derivatives of total energy with respect to electric fields
56 : : !! atom displacements,strain,...... all in cartesian coordinates
57 : : !! crystal<crystal_t>=Crystalline structure info.
58 : : !! iblok= blok number in DDB file
59 : : !! iout=out file number
60 : : !! mpert=maximum number of ipert
61 : : !! msize=Maximum size of dynamical matrices and other perturbations (ddk, dde...)
62 : : !! natom=number of atoms in unit cell
63 : : !! nblok=number of total bloks in DDB file
64 : : !! prt_internalstr=if 2 or higher, print force and displacement internal strain,
65 : : !! if 1, print only force internal strain,
66 : : !! if 0, do not print internal strain.
67 : : !!
68 : : !! OUTPUT
69 : : !! instrain=force response internal strain tensor
70 : : !!
71 : : !! NOTES
72 : : !! In output of internal strain tensor,column runs from strain1 to
73 : : !! strain6(in Voigt notation),row runs from atom1x,atom1y,atom1z,atom2x,.......
74 : : !! sum rule is applied on the internal strain tensor
75 : : !!
76 : : !! SOURCE
77 : :
78 : 35 : subroutine ddb_internalstr(asr,&
79 : : !&crystal,&
80 : 35 : & blkval,&
81 : : !&asrq0,&
82 : 35 : & d2asr,iblok,instrain,iout,mpert,&
83 : : !&msize,&
84 : : natom,nblok,prt_internalstr)
85 : :
86 : : !Arguments----------------------------------------------
87 : : !scalars
88 : : integer,intent(in) :: asr,iblok,iout,mpert,natom,nblok,prt_internalstr
89 : : !integer,intent(in) :: msize
90 : : !type(crystal_t),intent(in) :: crystal
91 : : !type(asrq0_t),intent(inout) :: asrq0
92 : : !arrays
93 : : real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
94 : : real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok)
95 : : real(dp),intent(out) :: instrain(3*natom,6)
96 : :
97 : : !Local variables------------------------------------
98 : : !scalars
99 : : integer :: idirA,idirB,ier,ii1,ipertA,ipertB,ivarA,ivarB
100 : : character(len=500) :: direction,message
101 : : !arrays
102 : 70 : real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom)
103 : 70 : real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2)
104 : 70 : real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3)
105 : 70 : real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom),deviation(3,6)
106 : 70 : real(dp) :: eigval(3*natom-3),eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3)
107 : 70 : real(dp) :: eigvecp(2,3*natom,3*natom),instrain_dis(6,3*natom)
108 : 70 : real(dp) :: kmatrix(3*natom,3*natom),zhpev1(2,2*3*natom-4)
109 : 70 : real(dp) :: zhpev1p(2,2*3*natom-1),zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2)
110 : 70 : real(dp) :: d2cart(2,3*natom,3*natom)
111 : :
112 : : !***************************************************************
113 : :
114 : : !extract internal strain from DDB matrix data
115 : :
116 [ + + ]: 203 : do ipertA=1,natom
117 [ + + ]: 707 : do idirA=1,3
118 : 504 : ivarA=idirA+3*(ipertA-1)
119 [ + + ]: 3696 : do ivarB=1,6
120 [ + + ]: 3024 : if(ivarB<=3) then
121 : 1512 : idirB=ivarB
122 : 1512 : ipertB=natom+3
123 : : ! for the diagonal modulus
124 : : else if(ivarB>3) then
125 : 1512 : idirB=ivarB-3
126 : 1512 : ipertB=natom+4
127 : : ! for the shear modulus
128 : : end if
129 : 3528 : instrain(ivarA,ivarB)=(-1.0_dp)*blkval(1,idirA,ipertA,idirB,ipertB,iblok)
130 : : ! write(std_out,'(es16.6)')blkval(1,idirA,ipertA,idirB,ipertB,iblok)
131 : : end do
132 : : end do
133 : : end do
134 : : !according to the definition there is a minus sign before the second derivative
135 : :
136 : : !apply sum rule to the internal strain tensor
137 : 35 : deviation(:,:)=zero
138 [ + + ]: 245 : do ivarB=1,6
139 [ + + ]: 3269 : do ivarA=1,3*natom
140 [ + + ]: 3024 : if(mod(ivarA,3)==0)then
141 : 1008 : deviation(1,ivarB)=deviation(1,ivarB)+instrain(ivarA,ivarB)
142 : : end if
143 [ + + ]: 3024 : if(mod(ivarA,3)==1)then
144 : 1008 : deviation(2,ivarB)=deviation(2,ivarB)+instrain(ivarA,ivarB)
145 : : end if
146 [ + + ]: 3234 : if(mod(ivarA,3)==2)then
147 : 1008 : deviation(3,ivarB)=deviation(3,ivarB)+instrain(ivarA,ivarB)
148 : : end if
149 : : end do
150 : : end do
151 : :
152 [ + + ]: 245 : do ivarB=1,6
153 [ + + ]: 3269 : do ivarA=1,3*natom
154 [ + + ]: 3024 : if(mod(ivarA,3)==0)then
155 : 1008 : instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(1,ivarB)/natom
156 : : end if
157 [ + + ]: 3024 : if(mod(ivarA,3)==1)then
158 : 1008 : instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(2,ivarB)/natom
159 : : end if
160 [ + + ]: 3234 : if(mod(ivarA,3)==2)then
161 : 1008 : instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(3,ivarB)/natom
162 : : end if
163 : : end do
164 : : end do
165 : : !ending the sum rule
166 : :
167 : : !print the force response internal strain constants into the output file
168 [ + - ]: 35 : if(prt_internalstr>0)then
169 : 35 : write(message,'(a,a,a,a)')ch10,&
170 : 70 : & ' Force-response internal strain tensor','(Unit:Hartree/bohr)',ch10
171 : 35 : call wrtout(std_out,message,'COLL')
172 : 35 : call wrtout(iout,message,'COLL')
173 : :
174 : 35 : write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
175 : 70 : & 'strainyy','strainzz','strainyz','strainxz','strainxy'
176 : 35 : call wrtout(std_out,message,'COLL')
177 [ + + ]: 539 : do ii1=1,3*natom
178 [ + + ]: 504 : if(mod(ii1,3)==1)then
179 : 168 : direction='x'
180 [ + + ]: 336 : elseif(mod(ii1,3)==2)then
181 : 168 : direction='y'
182 : : elseif(mod(ii1,3)==0)then
183 : 168 : direction='z'
184 : : end if
185 : 504 : write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,' ',direction,&
186 : 504 : & instrain(ii1,1),instrain(ii1,2),instrain(ii1,3),&
187 : 1008 : & instrain(ii1,4),instrain(ii1,5),instrain(ii1,6)
188 : 539 : call wrtout(std_out,message,'COLL')
189 : : end do
190 : : endif
191 : :
192 : : !now write into the ddb output file
193 : 35 : write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
194 : 70 : & 'strainyy','strainzz','strainyz','strainxz','strainxy'
195 : 35 : call wrtout(iout,message,'COLL')
196 [ + + ]: 539 : do ii1=1,3*natom
197 [ + + ]: 504 : if(mod(ii1,3)==1)then
198 : 168 : direction='x'
199 [ + + ]: 336 : elseif(mod(ii1,3)==2)then
200 : 168 : direction='y'
201 : : elseif(mod(ii1,3)==0)then
202 : 168 : direction='z'
203 : : end if
204 : 504 : write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,' ',direction,&
205 : 504 : & instrain(ii1,1),&
206 : 504 : & instrain(ii1,2),instrain(ii1,3),&
207 : 1008 : & instrain(ii1,4),instrain(ii1,5),instrain(ii1,6)
208 : 539 : call wrtout(iout,message,'COLL')
209 : : end do
210 : :
211 : : ! ----------------------------------------------------------------------------------------
212 : :
213 : : !Try to get the displacement response internal strain tensor
214 : : !first need the inverse of force constant matrix
215 [ + + + + : 10293 : d2cart(1,:,:) = RESHAPE(blkval(1,1:3,1:natom,1:3,1:natom,iblok), SHAPE=[3*natom,3*natom])
+ + ]
216 [ + + + + ]: 10223 : d2cart(2,:,:) = zero
217 : :
218 : : !Eventually impose the acoustic sum rule
219 : : !FIXME: this might depend on ifcflag: impose that it is 0 or generalize
220 : 35 : call asria_corr(asr,d2asr,d2cart,natom,natom)
221 : : !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
222 [ + + + + ]: 10223 : kmatrix = d2cart(1,:,:)
223 [ + + + + ]: 10223 : Apmatr(:,:)=kmatrix(:,:)
224 : :
225 : : !DEBUG
226 : : !write(std_out,'(/,a,/)')'the force constant matrix'
227 : : !do ivarA=1,3*natom
228 : : !write(std_out,'(/)')
229 : : !do ivarB=1,3*natom
230 : : !write(std_out,'(es16.6)')kmatrix(ivarB,ivarA)
231 : : !end do
232 : : !end do
233 : : !ENDDEBUG
234 : :
235 [ + + + + ]: 10223 : Nmatr(:,:)=0.0_dp
236 [ + + ]: 539 : do ivarA=1,3*natom
237 [ + + ]: 10223 : do ivarB=1,3*natom
238 [ + + + + ]: 9684 : if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0)then
239 : 1076 : Nmatr(ivarA,ivarB)=one
240 : : end if
241 [ + + + + ]: 9684 : if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1)then
242 : 1076 : Nmatr(ivarA,ivarB)=one
243 : : end if
244 [ + + + + ]: 10188 : if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2)then
245 : 1076 : Nmatr(ivarA,ivarB)=one
246 : : end if
247 : : end do
248 : : end do
249 : :
250 : : !DEBUG
251 : : !do ivarA=1,3*natom
252 : : !write(std_out,'(/)')
253 : : !do ivarB=1,3*natom
254 : : !write(std_out,'(es16.6)')Nmatr(ivarB,ivarA)
255 : : !end do
256 : : !end do
257 : : !ENDDEBUG
258 : :
259 [ + - ]: 35 : if (natom > 1) then
260 : :
261 : : !starting the pseudoinverting processes
262 : : !then get the eigenvectors of the big matrix,give values to matrixBp
263 [ + + + + ]: 15317 : Bpmatr=0.0_dp
264 : : ii1=1
265 [ + + ]: 539 : do ivarA=1,3*natom
266 [ + + ]: 5633 : do ivarB=1,ivarA
267 : 5094 : Bpmatr(1,ii1)=Nmatr(ivarB,ivarA)
268 : 5598 : ii1=ii1+1
269 : : end do
270 : : end do
271 : :
272 : : !Bpmatr(2,:) is the imaginary part of the force matrix
273 : : !then call the subroutines CHPEV and ZHPEV to get the eigenvectors
274 : 35 : call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier)
275 [ - + ]: 35 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
276 : :
277 : : !DEBUG
278 : : !the eigenval and eigenvec
279 : : !write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors'
280 : : !do ivarA=1,3*natom
281 : : !write(std_out,'(/)')
282 : : !write(std_out,'(es16.6)')eigvalp(ivarA)
283 : : !end do
284 : : !do ivarA=1,3*natom
285 : : !write(std_out,'(/)')
286 : : !do ivarB=1,3*natom
287 : : !write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA)
288 : : !end do
289 : : !end do
290 : : !ENDDEBUG
291 : :
292 : : !Then do the multiplication to get the reduced matrix,in two steps
293 : : !After this the force constant matrix is decouple in two bloks,
294 : : !acoustic and optical ones
295 [ + + + + ]: 10223 : Cpmatr(:,:)=0.0_dp
296 [ + + ]: 539 : do ivarA=1,3*natom
297 [ + + ]: 10223 : do ivarB=1,3*natom
298 [ + + ]: 311832 : do ii1=1,3*natom
299 : 311328 : Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*Apmatr(ii1,ivarB)
300 : : end do
301 : : end do
302 : : end do
303 : :
304 [ + + + + ]: 10223 : Apmatr(:,:)=0.0_dp
305 [ + + ]: 539 : do ivarA=1,3*natom
306 [ + + ]: 10223 : do ivarB=1,3*natom
307 [ + + ]: 311832 : do ii1=1,3*natom
308 : 311328 : Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*eigvecp(1,ii1,ivarB)
309 : : end do
310 : : end do
311 : : end do
312 : :
313 : : !DEBUG
314 : : !the blok diago
315 : : !write(std_out,'(/,a,/)')'matrixAp'
316 : : !do ivarA=1,3*natom
317 : : !write(std_out,'(/)')
318 : : !do ivarB=1,3*natom
319 : : !write(std_out,'(es16.6)')Apmatr(ivarA,ivarB)
320 : : !end do
321 : : !end do
322 : : !ENDDEBUG
323 : :
324 : : !Check the last three eigenvalues whether too large or not
325 : : ivarB=0
326 [ + + ]: 140 : do ivarA=3*natom-2,3*natom
327 [ - + ]: 140 : if (ABS(Apmatr(ivarA,ivarA))>tol6)then
328 : 0 : ivarB=1
329 : : end if
330 : : end do
331 : :
332 [ - + ]: 35 : if(ivarB==1)then
333 : 0 : write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,&
334 : 0 : & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,&
335 : 0 : & ' are too large at Gamma point.',ch10,&
336 : 0 : & ' Increase cutoff energy or k-points sampling.',ch10,&
337 : 0 : & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom)
338 : 0 : ABI_WARNING(message)
339 : 0 : call wrtout(iout,message,'COLL')
340 : : end if
341 : :
342 : : !Give the value of reduced matrix form Apmatr to Amatr
343 [ + + ]: 434 : do ivarA=1,3*natom-3
344 [ + + ]: 7409 : do ivarB=1,3*natom-3
345 : 7374 : Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)
346 : : end do
347 : : end do
348 : :
349 : : !Now the reduced matrix is in the matrixA, the convert it
350 : : !first give the give the value of matixB from matrixA
351 : : ii1=1
352 [ + + ]: 434 : do ivarA=1,3*natom-3
353 [ + + ]: 4121 : do ivarB=1,ivarA
354 : 3687 : Bmatr(1,ii1)=Amatr(ivarB,ivarA)
355 : 4086 : ii1=ii1+1
356 : : end do
357 : : end do
358 [ + + ]: 3722 : Bmatr(2,:)=0.0_dp
359 : :
360 : : !Call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues
361 : 35 : call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier)
362 [ - + ]: 35 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
363 : :
364 : : !Check the unstable phonon modes, if the first is negative then print
365 : : !warning message
366 [ + + ]: 35 : if(eigval(1)<-1.0*tol8)then
367 : 8 : write(message,'(9a)') ch10,&
368 : 8 : & ' --- !WARNING',ch10,&
369 : 8 : & ' Unstable eigenvalue detected in force constant matrix at Gamma point',ch10,&
370 : 8 : & ' The system under calculation is physically unstable.',ch10,&
371 : 16 : & ' ---',ch10
372 : 8 : call wrtout(std_out,message,'COLL')
373 : : end if
374 : :
375 : : !Do the matrix mutiplication to get pseudoinverse inverse matrix
376 [ + + + + ]: 7409 : Cmatr(:,:)=0.0_dp
377 [ + + + + ]: 7409 : Amatr(:,:)=0.0_dp
378 [ + + ]: 434 : do ivarA=1,3*natom-3
379 : 434 : Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA)
380 : : end do
381 : :
382 [ + + ]: 434 : do ivarA=1,3*natom-3
383 [ + + ]: 7409 : do ivarB=1,3*natom-3
384 [ + + ]: 234525 : do ii1=1,3*natom-3
385 : 234126 : Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*Cmatr(ii1,ivarB)
386 : : end do
387 : : end do
388 : : end do
389 : :
390 : :
391 : : !The second multiplication
392 [ + + + + ]: 7409 : Cmatr(:,:)=0.0_dp
393 [ + + ]: 434 : do ivarA=1,3*natom-3
394 [ + + ]: 7409 : do ivarB=1,3*natom-3
395 [ + + ]: 234525 : do ii1=1,3*natom-3
396 : 234126 : Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+ Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1)
397 : : end do
398 : : end do
399 : : end do
400 : :
401 : : !DEBUG
402 : : !write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix'
403 : : !do ivarA=1,3*natom
404 : : !write(std_out,'(/)')
405 : : !do ivarB=1,3*natom
406 : : !write(std_out,'(es16.6)')Cmatr(ivarA,ivarB)
407 : : !end do
408 : : !end do
409 : : !ENDDEBUG
410 : :
411 : : !So now the inverse of the reduced matrix is in the matrixC
412 : : !now do another mutilplication to get the pseudoinverse of the original
413 [ + + + + ]: 10223 : Cpmatr(:,:)=0.0_dp
414 [ + + + + ]: 10223 : Apmatr(:,:)=0.0_dp
415 [ + + ]: 434 : do ivarA=1,3*natom-3
416 [ + + ]: 7409 : do ivarB=1,3*natom-3
417 : 7374 : Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)
418 : : end do
419 : : end do
420 : :
421 : : !Now times the eigvecp
422 [ + + ]: 539 : do ivarA=1,3*natom
423 [ + + ]: 10223 : do ivarB=1,3*natom
424 [ + + ]: 311832 : do ii1=1,3*natom
425 : : Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*&
426 : 311328 : & Cpmatr(ii1,ivarB)
427 : : end do
428 : : end do
429 : : end do
430 [ + + + + ]: 10223 : Cpmatr(:,:)=0.0_dp
431 [ + + ]: 539 : do ivarA=1,3*natom
432 [ + + ]: 10223 : do ivarB=1,3*natom
433 [ + + ]: 311832 : do ii1=1,3*natom
434 : 311328 : Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+ Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1)
435 : : end do
436 : : end do
437 : : end do
438 : :
439 : : !Now the inverse is in Cpmatr
440 [ + + + + ]: 10223 : kmatrix(:,:)=Cpmatr(:,:)
441 : : !transfer the inverse of k-matrix back to the k matrix
442 : : !so now the inverse of k matrix is in the kmatrix
443 : : !ending the part for pseudoinversing the K matrix
444 : :
445 : : !Now do simple mulplication to obtain the displacement response
446 : : !internal strain tensor
447 [ + + + + ]: 3563 : instrain_dis(:,:)=0.0_dp
448 [ + + ]: 245 : do ivarA=1,6
449 [ + + ]: 3269 : do ivarB=1,3*natom
450 [ + + ]: 61338 : do ii1=1,3*natom
451 : : instrain_dis(ivarA,ivarB)=instrain_dis(ivarA,ivarB)+&
452 : 61128 : & instrain(ii1,ivarA)*kmatrix(ii1,ivarB)
453 : : end do
454 : : end do
455 : : end do
456 : :
457 : : else
458 [ # # # # ]: 0 : instrain_dis(:,:)=0.0_dp
459 : : end if
460 : :
461 : : !Print out the results
462 [ + + ]: 35 : if(prt_internalstr>1)then
463 : 6 : write(message,'(a,a,a,a)')ch10,&
464 : 12 : & ' Displacement-response internal strain ', 'tensor (Unit:Bohr)',ch10
465 : 6 : call wrtout(std_out,message,'COLL')
466 : 6 : call wrtout(iout,message,'COLL')
467 : 6 : write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
468 : 12 : & 'strainyy','strainzz','strainyz','strainxz','strainxy'
469 : 6 : call wrtout(std_out,message,'COLL')
470 : 6 : call wrtout(iout,message,'COLL')
471 [ + + ]: 54 : do ivarA=1,3*natom
472 [ + + ]: 48 : if(mod(ivarA,3)==1)then
473 : 16 : direction='x'
474 [ + + ]: 32 : elseif(mod(ivarA,3)==2)then
475 : 16 : direction='y'
476 : : elseif(mod(ivarA,3)==0)then
477 : 16 : direction='z'
478 : : end if
479 : 48 : write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ivarA-1)/3)+1,' ',direction,&
480 : 48 : & instrain_dis(1,ivarA),instrain_dis(2,ivarA),&
481 : 48 : & instrain_dis(3,ivarA),instrain_dis(4,ivarA),instrain_dis(5,ivarA),&
482 : 96 : & instrain_dis(6,ivarA)
483 : 48 : call wrtout(std_out,message,'COLL')
484 : 54 : call wrtout(iout,message,'COLL')
485 : : end do
486 : : endif
487 : :
488 : 35 : end subroutine ddb_internalstr
489 : : !!***
490 : :
491 : : end module m_ddb_internalstr
492 : : !!***
|
