-
Notifications
You must be signed in to change notification settings - Fork 6
gulp af
###links http://fortran90.org/src/best-practices.html#modules-and-programs
http://www.csee.umbc.edu/~squire/fortranclass/summary.shtml#Stru
#QUESTIONS: 01/06/12 ##Q
! use constants,
use current
use general, only : bfactor
use iochannels
use parallel
implicit none
###question
where are the modules at?
####answer
Refer to file:
~/Src/modules.f90
##Q
subroutine thermalconductivity(mcv,derv2,eigr,Sij,freq,nphonatc,ncfoc,nphonatptr,maxd2)
###question why is Sij listed with the subroutine definition, when it seems to me to be more of a "local" variable since its contents are defined within the thermalconductivtiy.f90? ####answer It could be local as things stand. However Sij is the result of the subroutine and so it is conceivable that it might be needed as a return argument. For example, if fitting of Sij was allowed then this might be the case. There is also the possibility that memory overlays can be used (i.e. an array that is used elsewhere in the calling routine "phonon" could be passed as the Sij workspace). ##Q
!
! Passed variables
!
integer(i4), intent(in) :: nphonatptr(*)
integer(i4), intent(in) :: maxd2
integer(i4), intent(in) :: mcv
integer(i4), intent(in) :: ncfoc
integer(i4), intent(in) :: nphonatc
real(dp), intent(inout) :: derv2(maxd2,*)
real(dp), intent(in) :: eigr(maxd2,*)
real(dp), intent(out) :: Sij(mcv,*)
real(dp), intent(in) :: freq(*)
###question where do these variables get passed in from the main part of the code? this is one of the challenges for me, tracing functionality back to the part of the code which reads in the gulp.gin. ####answer These get passed in from "phonon" as the calling master phonon calculation routine. The best thing is rather than have you track things back in the code from the input file that I tell you what the variables are in a physical sense:
maxd2 - this is the left-hand dimension of the arrays derv2 and eigr
mcv - this is the left-hand dimension of Sij and is 3 x ncfoc
ncfoc - this is the number of fully occupied cores in the system. In the simplest case this is just equal to the number of atoms in the unit cell (numat). In cases where there is a shell model being used or there are partially occupied sites then it can be < numat.
nphonatc - this is the number of cores in the system. Again, in the simplest case it's just equal to numat.
nphonatptr(*) - this is an array of integers that point from the number of the core (from 1 to nphonatc) - to the number of full coordinate list site (that includes shells). For the simple case, then nphonatptr(i) = i.
derv2(maxd2,) - this is the dynamical matrix (i.e. mass weighted second derivatives) as a 3ncfoc x 3*ncfoc matrix
eigr(maxd2,) - these are the eigenvectors as a 3ncfoc x 3*ncfoc matrix
freq() - this is a linear array containing the frequencies in ascending order of the 3ncfoc modes
###question what is the meaning of the syntax "derv2(maxd2,*)"? what will the size and shape of the array be? ####answer
I think I answered my own question, line 164:
!
! Copy results back to Sij
!
Sij(1:mcv,1:mcv) = derv2(1:mcv,1:mcv)
what is the variable "mcv"? It must be the number of modes in the system, 3num_atoms_total. Based on this, how do I know to label the passed array as "derv2(maxd2,)"? ##Q mcv is indeed 3 x num_atoms_total (numat in GULP). Not sure quite what the second question is asking?
!
! Allocate local array ldone to avoid duplicate multiplies in case of partial occupancy
!
ncfoc2 = ncfoc*(ncfoc+1)/2
allocate(ldone(ncfoc2),stat=status)
if (status/=0) call outofmemory('thermalconductivity','ldone')
ldone(1:ncfoc2) = .false.
###question what does "partial occupancy" mean? I guess ldone(ncfoc2) is an array of logical ".true." or ".false." of length ncfoc2, based on whether the array was allocated correctly? ####answer Because 2 or more atoms in the full cell can be linked as being part of the same fully occupied site (due to partial occupancies) then their matrix elements have to be combined. The flag ldone indicates whether this site has been addressed already or not. ##Q
!
! Scale dynamical matrix elements by minimum image nearest distance between sites
!
do i = 1,nphonatc
ii = nphonatptr(i)
ix = 3*ii - 2
iy = ix + 1
iz = ix + 2
do j = 1,i-1
jj = nphonatptr(j)
ind = ii*(ii-1)/2 + jj
if (.not.ldone(ind)) then
jx = 3*jj - 2
jy = jx + 1
jz = jx + 2
###question what is the meaning of the variable "nphonatc" and "nphonatptr"? I guess it they are related to either the number of modes, or number of atoms, or both. ####answer This is covered in an earlier answer. ##Q
!
! Compute initial vector
!
xd = xclat(jj) - xclat(ii)
yd = yclat(jj) - yclat(ii)
zd = zclat(jj) - zclat(ii)
###question the variables "{x,y,z}clat" are defined as local variables, but don't seem to be allocated or initialized anywhere before these lines? Does this mean all values are = 0? ####answer {x,y,z}clat are actually global variables that come from modules. They are the x, y and z, respectively, Cartesian coordinates of the atoms in Angstroms.
#NEW
#QUESTIONS: 01/08/12
##Q
!
! Create inverse frequency factors while trapping translations and imaginary modes
!
allocate(freqinv(mcv),stat=status)
if (status/=0) call outofmemory('thermalconductivity','freqinv')
nfreqmin = 0
do i = 1,mcv
if (freq(i).gt.1.0_dp) then
if (nfreqmin.eq.0) nfreqmin = i
freqinv(i) = 1.0_dp/sqrt(2.0_dp*freq(i))
else
freqinv(i) = 1.0_dp
endif
enddo
###question do you store imag freq as negative numbers also internally in the code? I know neg freq == imag in the output.gout. Also, why do you check "freq(i).gt.1.0_dp" and not "freq(i).gt.0.0_dp"? ####answer GULP does indeed store imaginary frequencies as a negative number for convenience. The reason for the check that the frequency > 1 is to eliminate the divide by zero error that would occur for the translational modes in the lines that follow. This might be a bit crude, but is unlikely that most systems would have a genuine frequency that was below 1 cm^-1. ##Q ###question What are the units of "freq" internally?
####answer The units of freq are in wave numbers (cm^-1) on return from pdiag(g). ##Q ###question could you add an option to output the eigenvectors like the freq:
output freq freq.gout
output eig eig.gout ####answer If you add the keyword "eigen" to the input file then you should get a print out of the eigenvectors in the .gout file. If you want them in a separate file then I could look to add this as an option, but would obviously take longer.
##Q ###question And also an option to output the D_i:
output Di Di.gout ####answer The Di values are given in the output of the program. However, I can again add an option to write to a separate file if you let me know what would be the preferred format.
##Q
###question
####answer
#result
##GULP 4.0.5
note the units are freq[cm^{-1}] and Dij[not sure yet]
##GULP 4.0.7-units
###thermalconductivity_gale_022413
source code:
###thermalconductivity
I changed your defn of the lorentzian, which seems to be mising a factor of pi, to:
dwij = (1.0/pi)*(bfactor*dwavg)/( (freq(j) - freq(i))**2 + (bfactor*dwavg)**2 )
and this result looks identical to what my matlab code gets, except for the unit conversion.
https://raw.github.com/jasonlarkin/disorder/master/matlab/gulp/lj_amor_4.0.7_units-lor.jpg
source code:
here is the input.gin for that structure if you want to play with it to figure out the units:
##correct
note the units are LJ units for both
###units obviously the unit conversion above is no good, here are my factors:
function Di = m_gulp_af_lj_readDi(path,name)
con = m_constant; lj = m_lj;
Di = load([path name]);
%gets to m^2/s
factor.gulp = con.eV2J^2*con.avog^2/con.ang2m^2/con.c
factor =
gulp: 3.105287902711237e+19
con.eV2J
ans =
1.602176460000000e-19
con.avog
ans =
6.022141290000000e+23
con.ang2m
ans =
1.000000000000000e-10
con.c
ans =
2.997924580000019e+10
%gets to lj units
factor.lj = lj.tau/lj.sigma^2
Di_convert(:,3) = Di(:,3)*factor.gulp*factor.lj;
Di_convert(:,2) = Di(:,2)*con.c*lj.tau;
loglog(Di_convert(:,2),Di_convert(:,3),'.')
end
##GULP 4.0.7-lor
I changed the following:
real(dp) :: pi
pi = 3.14159265359
!
constant = 1.0d10*pi/12.0_dp ! 1/3 convoluted with 1/2 squared from A7
!
! Scale by constants and inverse frequency squared - factor of third is for averaging over directions
!
Di(i) = Di(i) + Di_loc*constant/(freq(i)**2)
This is alot closer:
#conv
The updated version of the thermalconductivity.f90:
https://github.com/jasonlarkin/disorder/blob/master/matlab/gulp/thermalconductivity.f90
I changed the broadening to the average freq spacing, and now the variable broaden controls what multiple of dwavg the Lor is broadened by:
! DEBUG
!
! Find mean level spacing
!
dwavg = 0.0_dp;
do i = nfreqmin,mcv-1
dwavg = freq(i+1) - freq(i) + dwavg
enddo
dwavg = dwavg/(mcv-1 - nfreqmin)
! DEBUG
!
! Sum over coupling with mode j weighted by Lorentzian factor
!
do j = nfreqmin,i-1
dwij = (bfactor*dwavg)/(1 + ((bfactor*dwavg)*(freq(j) - freq(i)))**2)
Di = Di + dwij*Sij(j,i)**2
enddo
do j = i+1,mcv
dwij = (bfactor*dwavg)/(1 + ((bfactor*dwavg)*(freq(j) - freq(i)))**2)
Di = Di + dwij*Sij(j,i)**2
enddo
I also added a print out of the values of i,j,rij:
! DEBUG
!
! Write out a file with rij
!
open(91,file='rij.gout',status='unknown',form='formatted')
open(91,file='rij.gout',status='unknown',form='formatted')
write(91,'(i6,2x,i6,2x,f12.4)') &
&i,j,rij
! DEBUG
Here is a comparison for the following structure:
4 4 20.000 20.000 20.000
1 1 0.0000 0.0000 0.0000
2 1 0.7818 0.7818 0.0000
3 1 0.7818 0.0000 0.7818
4 1 0.0000 0.7818 0.7818
##GULP
The results are in this folder:
https://github.com/jasonlarkin/disorder/tree/master/matlab/gulp/conv
##matlab
Here is the matlab code used to produce the below:
https://github.com/jasonlarkin/disorder/blob/master/matlab/m_af_lj.m
###freq
0+2.0197i
0+2.0197i
0+2.0197i
0+7.6649e-08i
1.0875e-07
2.9761e-07
8.6421
8.6421
12.387
12.387
12.387
17.635
####gulp
-5.003599416318
-5.003599416318
-5.003599416318
-0.000000401238
-0.000000215338
0.000000055440
21.410294901138
21.410294901138
30.689370515317
30.689370515317
30.689370515317
43.688796417065
which converts to lj units as:
-2.0195
-2.0195
-2.0195
-0.0000
-0.0000
0.0000
8.6414
8.6414
12.3866
12.3866
12.3866
17.6333
###Di
0+2.0197i -0
0+2.0197i -0
0+2.0197i -0
0+7.6649e-08i -9.5417e-16-6.604e-30i
1.0875e-07 1.5531e-15-3.7303e-30i
2.9761e-07 1.1245e-16-3.4781e-31i
8.6421 0.052754
8.6421 0.052754
12.387 0.092533
12.387 0.092533
12.387 0.092533
17.635 0.11164
###eigvec
For the eigenvectors, I suppose there is degeneracy (or near degeneracy), so the eigenvectors don't match up exactly:
####matlab
0.45418 0.20003 0.060894 0.47816 -0.0037812 -0.14611 -0.20412 0.35355 0.28183 0.35721 0.2073 -0.28868
-0.36588 0.26938 -0.20873 0.14424 -0.068707 0.47379 -0.20412 -0.35355 0.28336 0.2939 -0.28867 -0.28868
-0.088302 -0.46941 0.14783 -0.023661 -0.49524 -0.064615 0.40825 -3.4235e-16 -0.12093 0.4812 -0.061803 -0.28868
-0.20219 0.38301 0.24985 0.47816 -0.0037812 -0.14611 0.20412 -0.35355 -0.40429 0.18731 0.22687 0.28868
0.29049 0.086404 -0.39768 0.14424 -0.068707 0.47379 0.20412 0.35355 -0.40275 0.12399 -0.2691 0.28868
0.088302 0.46941 -0.14783 -0.023661 -0.49524 -0.064615 0.40825 -6.1542e-16 0.12093 -0.4812 0.061803 -0.28868
0.20219 -0.38301 -0.24985 0.47816 -0.0037812 -0.14611 0.20412 -0.35355 0.40429 -0.18731 -0.22687 0.28868
0.36588 -0.26938 0.20873 0.14424 -0.068707 0.47379 -0.20412 -0.35355 -0.28336 -0.2939 0.28867 -0.28868
0.16369 0.11363 0.45858 -0.023661 -0.49524 -0.064615 -0.40825 5.9968e-16 0.0015339 -0.063314 -0.49597 0.28868
-0.45418 -0.20003 -0.060894 0.47816 -0.0037812 -0.14611 -0.20412 0.35355 -0.28183 -0.35721 -0.2073 -0.28868
-0.29049 -0.086404 0.39768 0.14424 -0.068707 0.47379 0.20412 0.35355 0.40275 -0.12399 0.2691 0.28868
-0.16369 -0.11363 -0.45858 -0.023661 -0.49524 -0.064615 -0.40825 1.5558e-16 -0.0015339 0.063314 0.49597 0.28868
####GULP
0.080147 -0.268189 0.188042 0.188701 -0.376743 -0.188042
-0.188701 0.268189 0.456890 -0.080147 0.376743 -0.456890
-0.036137 0.380558 -0.344420 0.463013 -0.118593 0.344420
-0.463013 -0.380558 0.082456 0.036137 0.118593 -0.082456
-0.492210 0.182347 0.309862 -0.003267 -0.306595 -0.309862
0.003267 -0.182347 -0.185615 0.492210 0.306595 0.185615
-0.029092 0.033400 -0.498034 -0.029092 0.033400 -0.498034
-0.029092 0.033400 -0.498034 -0.029092 0.033400 -0.498034
0.447081 -0.220108 -0.040877 0.447081 -0.220108 -0.040877
0.447081 -0.220108 -0.040877 0.447081 -0.220108 -0.040877
0.221973 0.447702 0.017058 0.221973 0.447702 0.017058
0.221973 0.447702 0.017058 0.221973 0.447702 0.017058
-0.353553 0.353553 -0.000000 0.353553 -0.353553 -0.000000
0.353553 0.353553 0.000000 -0.353553 -0.353553 0.000000
-0.204124 -0.204124 0.408248 0.204124 0.204124 0.408248
0.204124 -0.204124 -0.408248 -0.204124 0.204124 -0.408248
0.231515 0.451818 0.402551 -0.049267 0.171036 -0.402551
0.049267 -0.451818 0.220303 -0.231515 -0.171036 -0.220303
0.363723 -0.081714 0.223299 0.305013 -0.140424 -0.223299
-0.305013 0.081714 -0.445437 -0.363723 0.140424 0.445437
-0.253193 -0.197948 0.195168 0.393115 0.448361 -0.195168
-0.393115 0.197948 0.055246 0.253193 -0.448361 -0.055246
-0.288675 -0.288675 -0.288675 0.288675 0.288675 -0.288675
0.288675 -0.288675 0.288675 -0.288675 0.288675 0.288675
The below rij{x,y,z} and Sij{x,y,z} may be helpful
###LD.Dynam
before LD.Dynam = (1/2)*(LD.Dynam' + LD.Dynam)
75.7050 39.3821 39.3821 -38.3624 -39.3821 0 -38.3624 0 -39.3821 1.0198 0 0
39.3821 75.7050 39.3821 -39.3821 -38.3624 0 0 1.0198 0 0 -38.3624 -39.3821
39.3821 39.3821 75.7050 0 0 1.0198 -39.3821 0 -38.3624 0 -39.3821 -38.3624
-38.3624 -39.3821 0 75.7050 39.3821 -39.3821 1.0198 0 0 -38.3624 0 39.3821
-39.3821 -38.3624 0 39.3821 75.7050 -39.3821 0 -38.3624 39.3821 0 1.0198 0
0 0 1.0198 -39.3821 -39.3821 75.7050 0 39.3821 -38.3624 39.3821 0 -38.3624
-38.3624 0 -39.3821 1.0198 0 0 75.7050 -39.3821 39.3821 -38.3624 39.3821 0
0 1.0198 0 0 -38.3624 39.3821 -39.3821 75.7050 -39.3821 39.3821 -38.3624 0
-39.3821 0 -38.3624 0 39.3821 -38.3624 39.3821 -39.3821 75.7050 0 0 1.0198
1.0198 0 0 -38.3624 0 39.3821 -38.3624 39.3821 0 75.7050 -39.3821 -39.3821
0 -38.3624 -39.3821 0 1.0198 0 39.3821 -38.3624 0 -39.3821 75.7050 39.3821
0 -39.3821 -38.3624 39.3821 0 -38.3624 0 0 1.0198 -39.3821 39.3821 75.7050
after LD.Dynam = (1/2)*(LD.Dynam' + LD.Dynam)
75.7050 39.3821 39.3821 -38.3624 -39.3821 0 -38.3624 0 -39.3821 1.0198 0 0
39.3821 75.7050 39.3821 -39.3821 -38.3624 0 0 1.0198 0 0 -38.3624 -39.3821
39.3821 39.3821 75.7050 0 0 1.0198 -39.3821 0 -38.3624 0 -39.3821 -38.3624
-38.3624 -39.3821 0 75.7050 39.3821 -39.3821 1.0198 0 0 -38.3624 0 39.3821
-39.3821 -38.3624 0 39.3821 75.7050 -39.3821 0 -38.3624 39.3821 0 1.0198 0
0 0 1.0198 -39.3821 -39.3821 75.7050 0 39.3821 -38.3624 39.3821 0 -38.3624
-38.3624 0 -39.3821 1.0198 0 0 75.7050 -39.3821 39.3821 -38.3624 39.3821 0
0 1.0198 0 0 -38.3624 39.3821 -39.3821 75.7050 -39.3821 39.3821 -38.3624 0
-39.3821 0 -38.3624 0 39.3821 -38.3624 39.3821 -39.3821 75.7050 0 0 1.0198
1.0198 0 0 -38.3624 0 39.3821 -38.3624 39.3821 0 75.7050 -39.3821 -39.3821
0 -38.3624 -39.3821 0 1.0198 0 39.3821 -38.3624 0 -39.3821 75.7050 39.3821
0 -39.3821 -38.3624 39.3821 0 -38.3624 0 0 1.0198 -39.3821 39.3821 75.7050
###Sij #####rijx
0 0.7818 0.7818 0
-0.7818 0 0 -0.7818
-0.7818 0 0 -0.7818
0 0.7818 0.7818 0
#####Rijx
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0
#####vijx
0 0.0000 -0.0000 -0.2829 0.5418 -0.9398 2.8429 1.6414 5.6971 24.6446 17.2445 0.0000
-0.0000 -0.0000 -0.0000 0.1282 0.8807 0.3966 -13.4162 -7.7458 26.0129 -0.0352 -5.5945 0.0000
0.0000 0.0000 0 0.1104 0.4494 0.3495 22.8669 13.2022 14.5536 -3.0846 -5.4262 -0.0000
0.2829 -0.1282 -0.1104 0.0000 0.0000 0.0000 -11.3980 19.7419 12.2197 1.3943 -1.5710 -67.1178
-0.5418 -0.8807 -0.4494 -0.0000 -0.0000 -0.0000 0.0901 -0.1561 1.6203 -33.7517 -25.6339 0.5308
0.9398 -0.3966 -0.3495 -0.0000 0.0000 0.0000 3.4829 -6.0326 39.9471 5.4363 -4.4777 20.5096
-2.8429 13.4162 -22.8669 11.3980 -0.0901 -3.4829 -0.0000 -0.0000 -1.5007 5.2527 -7.0110 -0.0000
-1.6414 7.7458 -13.2022 -19.7419 0.1561 6.0326 0.0000 0.0000 2.5993 -9.0979 12.1434 -0.0000
-5.6971 -26.0129 -14.5536 -12.2197 -1.6203 -39.9471 1.5007 -2.5993 0.0000 0.0000 0.0000 8.4893
-24.6446 0.0352 3.0846 -1.3943 33.7517 -5.4363 -5.2527 9.0979 -0.0000 0.0000 0.0000 -29.7137
-17.2445 5.5945 5.4262 1.5710 25.6339 4.4777 7.0110 -12.1434 -0.0000 -0.0000 -0.0000 39.6601
-0.0000 -0.0000 0.0000 67.1178 -0.5308 -20.5096 0.0000 0.0000 -8.4893 29.7137 -39.6601 -0.0000
#####Sijx
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0.0000 + 0.0000i 0.0000 + 0.0000i 0 0 0 0 0 0
0 0 0 0.0000 - 0.0000i 0 0.0000 0 0 0 0 0 0
0 0 0 0.0000 - 0.0000i 0.0000 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0.0000 0.0179 0.2190 0.3902 0
0 0 0 0 0 0 0.0000 0 0.0536 0.6570 1.1705 0
0 0 0 0 0 0 0.0179 0.0536 0 0.0000 0.0000 0.3150
0 0 0 0 0 0 0.2190 0.6570 0.0000 0 0.0000 3.8597
0 0 0 0 0 0 0.3902 1.1705 0.0000 0.0000 0 6.8762
0 0 0 0 0 0 0 0 0.3150 3.8597 6.8762 0
#####rijy
0 0.7818 0 0.7818
-0.7818 0 -0.7818 0
0 0.7818 0 0.7818
-0.7818 0 -0.7818 0
#####Rijy
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0
#####vijy
0.0000 0.0000 -0.0000 0.9952 -0.1270 -0.3214 -9.5024 5.4862 1.6403 -17.6789 22.6348 0.0000
-0.0000 -0.0000 -0.0000 -0.2522 0.5641 0.1586 -21.9857 12.6934 13.6864 8.4309 -6.7134 -0.0000
0.0000 0.0000 -0.0000 -0.3339 -0.9554 -0.0369 -11.7178 6.7653 -27.0096 -1.4822 -5.7593 0.0000
-0.9952 0.2522 0.3339 -0.0000 0.0000 0.0000 -3.4382 -5.9551 39.6965 8.5602 0.4608 -20.2458
0.1270 -0.5641 0.9554 -0.0000 0 -0.0000 1.6378 2.8367 6.7822 -24.9050 33.1478 9.6442
0.3214 -0.1586 0.0369 0.0000 0.0000 0.0000 -11.2938 -19.5614 -11.1012 -6.2176 4.6667 -66.5042
9.5024 21.9857 11.7178 3.4382 -1.6378 11.2938 0.0000 0.0000 -1.5393 6.8443 5.4573 0.0000
-5.4862 -12.6934 -6.7653 5.9551 -2.8367 19.5614 -0.0000 -0.0000 -2.6661 11.8547 9.4524 -0.0000
-1.6403 -13.6864 27.0096 -39.6965 -6.7822 11.1012 1.5393 2.6661 -0.0000 -0.0000 0.0000 8.7074
17.6789 -8.4309 1.4822 -8.5602 24.9050 6.2176 -6.8443 -11.8547 -0.0000 0.0000 -0.0000 -38.7174
-22.6348 6.7134 5.7593 -0.4608 -33.1478 -4.6667 -5.4573 -9.4524 -0.0000 0.0000 0.0000 -30.8713
-0.0000 0.0000 -0.0000 20.2458 -9.6442 66.5042 -0.0000 0.0000 -8.7074 38.7174 30.8713 -0.0000
#####Sijy
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0.0000 + 0.0000i 0.0000 + 0.0000i 0 0 0 0 0 0
0 0 0 0.0000 - 0.0000i 0 0.0000 0 0 0 0 0 0
0 0 0 0.0000 - 0.0000i 0.0000 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0.0000 0.0188 0.3718 0.2364 0
0 0 0 0 0 0 0.0000 0 0.0564 1.1155 0.7092 0
0 0 0 0 0 0 0.0188 0.0564 0 0.0000 0.0000 0.3314
0 0 0 0 0 0 0.3718 1.1155 0.0000 0 0.0000 6.5532
0 0 0 0 0 0 0.2364 0.7092 0.0000 0.0000 0 4.1663
0 0 0 0 0 0 0 0 0.3314 6.5532 4.1663 0
#####rijz
0 0 0.7818 0.7818
0 0 0.7818 0.7818
-0.7818 -0.7818 0 0
-0.7818 -0.7818 0 0
#####Rijz
0 0 0 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000
0 0 0 0 0 0 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000 0.781800000000000
-0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0
-0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 -0.781800000000000 0 0 0 0 0 0
#####vijz
-0.0000 -0.0000 -0.0000 0.3496 0.0135 -0.2313 0.0000 28.5795 1.2841 -3.9566 10.6708 -0.0000
0.0000 0.0000 -0.0000 1.0527 -0.0850 0.2659 -0.0000 -7.9672 6.9694 -26.9328 10.5135 0.0000
0.0000 0.0000 0.0000 0.1949 0.1291 -1.0610 0.0000 -8.2276 -2.2883 12.3368 26.8856 -0.0000
-0.3496 -1.0527 -0.1949 0.0000 0 0.0000 -1.1280 0.0000 -9.0906 38.4663 15.2540 3.3212
-0.0135 0.0850 -0.1291 0.0000 0.0000 -0.0000 -23.6104 0.0000 1.0397 -3.6915 4.4006 69.5157
0.2313 -0.2659 1.0610 -0.0000 0.0000 0.0000 -3.0805 -0.0000 -4.6395 14.2077 -39.3138 9.0698
-0.0000 0.0000 -0.0000 1.1280 23.6104 3.0805 0.0000 0.0000 -17.2482 -4.2713 0.4919 0.0000
-28.5795 7.9672 8.2276 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
-1.2841 -6.9694 2.2883 9.0906 -1.0397 4.6395 17.2482 -0.0000 0.0000 0.0000 -0.0000 -48.7853
3.9566 26.9328 -12.3368 -38.4663 3.6915 -14.2077 4.2713 -0.0000 -0.0000 -0.0000 -0.0000 -12.0810
-10.6708 -10.5135 -26.8856 -15.2540 -4.4006 39.3138 -0.4919 -0.0000 0.0000 0.0000 -0.0000 1.3913
0.0000 -0.0000 0.0000 -3.3212 -69.5157 -9.0698 -0.0000 -0.0000 48.7853 12.0810 -1.3913 -0.0000
#####Sijz
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0.0000 + 0.0000i 0 0 0 0 0 0
0 0 0 0.0000 - 0.0000i 0 0.0000 0 0 0 0 0 0
0 0 0 0.0000 - 0.0000i 0.0000 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0.0000 2.3615 0.1448 0.0019 0
0 0 0 0 0 0 0.0000 0 0.0000 0.0000 0.0000 0
0 0 0 0 0 0 2.3615 0.0000 0 0.0000 0.0000 10.4044
0 0 0 0 0 0 0.1448 0.0000 0.0000 0 0.0000 0.6380
0 0 0 0 0 0 0.0019 0.0000 0.0000 0.0000 0 0.0085
0 0 0 0 0 0 0 0 10.4044 0.6380 0.0085 0