Re: [SIESTA-L] reg- PDOS

[apostnik]

Dear Amar -
sorry, the script is not sophisticated enough to pick up
the atom numbers from the list. (You are welcome to make changes...)
In your case - if there are many other N atoms which you don't need -
I'd suggest that you run the script separately for each atom number -
and then sum up the data from three files...

Best regards

Andrei Postnikov



> Dear Siesta Users,
>I want to plot PDOS for three Nitrogen atoms
> simultaneously (atom number- 15, 20 and 23) in zigzag boron nitride
> nanoribbon using fmpdos utility in SIESTA. Now what do I write when asked
> as under:
>
>
> amar@amar-Inspiron-660s:~/ZBN-48/zbnr-bands$ fmpdos
>   Input file name (PDOS):
> zbnr.PDOS
>   Output file name :
> 3-N.dat
>   Extract data for atom index (enter atom NUMBER, or 0 to select all),
>   or for all atoms of given species (enter its chemical LABEL):
> ?
> (15 20 23 or 15+20+23)
> ??
>
> regards ...
>
> Amar Bahadur
> KNIPSS, Sultanpur- 228118 (Uttar Pradesh) INDIA
> Mob- +91- 9451431428
>

Re: [SIESTA-L] Help regarding structural relaxations for AMoS2NR

[apostnik]

> Dear all,
>
> I am trying to relax an armchair MoS2 nanoribbon. My initial
> coordinates are taken from the MoS2 sheet. My .fdf file together with
> my intial positions are attached. The problem is that after about 70
> steps of CG relaxation, all the atoms are coming together and the
> distances between atoms are much shorter than the initial ones. The
> relaxed structure after about 70 steps is attached in elec.XV file. I
> have successfully done relaxations for the MoS2 sheet but it seems
> that I am not able to relax the nanoribbon. Any helps would be
> apprecited.
>
> Bests,
> Mohammad,
>

Dear Mohammad,
some suggestions:

1. If you start a complicated system and a mess comes out,
it might be a good idea to try to sort the things out. You have a bit
too much of everything: many atoms, 70 CG steps, empty supercell whose
dimensions you relax (that does not make sense anyway) etc.
How about trying just bulk MoS2 first?
Then, if it not not quite to expectations, you can pull different
levers: basis, cutoffs, XC etc. If everything works fine, you can
proceed to your nanoribbons or whatever.

2. If you expect help on a quite vague issue,
it might be a good idea to provide at least the output file.
There are not many cases when the error is seen right away from
the input file, and to expect the colleagues to repeat your calculations
just out of curiosity might be too optimistic.

Have a nice day,

Andrei Postnikov

Re: [SIESTA-L] Help with DM.InitSpin

[apostnik]

Dear Parsa Ravan,

I don't have experience with TranSiesta but it seems that
your questions concern plain Siesta definitions, so I'll try.


> Dear TranSiesta users,
> I am trying to do transport calculations for a
> ferromagnetic/molecule/ferromagnetic junction. I have three questions
> about
> the job and I would be very grateful if someone could help me:
>
> 1.By simply setting "SpinPolarized.true." the code runs very
> nicely
> but
>when I initialize the spin using the block:
>
>%block DM.InitSpin
>1   -
>%endblock DM.InitSpin
>
>the scf loop DOESN'T converge. I need to add the DM.InitSpin block
> because for the antiparralel configuration of the electrodes it is
> required
> that magnetic moment direction of the right electrode be reversed.

I don't know what is your whole calculation in this example but
please mind that the DM.InitSpin addresses atom by atom
and not all atoms of given species. So your above definition says
set maximal spin-down on atom 1 and leave the rest (if any)
non-magnetic. If this is indeed what you want it should work.
Otherwise, having already the spin-up case converged,
you can invert the spin density as you wish, using my tool
http://www.home.uni-osnabrueck.de/apostnik/Software/DMtune.tar.gz
(which may require counting basis functions in tricky cases).


> 2.What are the numbers in the 4th column of kgrid block?

They are displacements of the k-points...
It is a long story; you can check, e.g.,
Moreno and Soler  http://dx.doi.org/10.1103/PhysRevB.45.13891

>%block kgrid_Monkhorst_Pack
>   1   0   0  0.0
>   0   1   0  0.0
>   0   0  10  0.5
>%endblock Kgrid_Monkhorst_Pack
>The userguide says that they usually must be set 0.0 or 0.5. Why?

As unshifted mesh explicitly contains Gamma and other symmetric points,
whereas the shifted mesh misses the band extremities, but it contains
less inequivalent k-points than in unshifted mesh, at the same mesh density.
Why 0.0 and 0.5 ? - the reduced k-points ought to be symmetric somehow,
in order to be compact. Siesta does not use much of symmetry of k-points,
but at least time reversal symmetry allows to half the total number
of k-points. With an arbitrary displacement (not 0.0 or 0.5)
this won't be possible.


> 3.For the antiparallel alignment of the electrodes I need to
> differentiate
> between those atoms of the right and left electrodes that are included
> in the
> scattering region. I tried many forms of chemical species
> declarations,
> for
> example this one:
>
>NumberOfSpecies 5
>%block Chemical_Species_Label
>1   26  Fe1
>2   16  S
>36  C
>41  H
>5   26  Fe2
>%endblock Chemical_Species_Label
>
>but the code does not accept the syntax.

Sure; what is expected is just chemical label (I think).
It makes sense to differentiate chemically similar atoms only if you
attribute to them different pseudopots or different bases. Otherwise
different spin configurations etc. are managed addressing
the individual atom index, not the species index.

Best regards

Andrei Postnikov

Re: [SIESTA-L] Question regarding kpoints and their displacement

[apostnik]

Dear Mohammad -
well, this is a nice simple exercise in direct/reciprocal cells.

Let "a" be the C-C distance (as assumed in the the small-supercell input),
then the supercell translation vectors (in the graphene plane) are:
AA = [3a,0]; BB = [0, 2a*sqrt(3)]
and the reciprocal lattice vectors :
AA_inv = [2*pi/(3*a), 0]; BB_inv = [0, pi/(a*sqrt(3))].
Whereas translation vectors of the graphene primitive cell
in the same setting, not moving any atoms and not rotating anything,
are, for example:
A = a/2[3, sqrt(3)]; B=a/2[-3, sqrt(3)].
Their respective reciprocal vectors are:
A_inv = (2*pi)/(3*a)[1, sqrt(3)]; B_inv = (2*pi)/(3*a)[-1, sqrt(3)].
With these, you can draw the Brillouin zone
and check where the K points are.
There are six around Gamma:
(2*pi)/(3*a)[(+/-)1, (+/-)1/sqrt(3)]
and
(2*pi)/(3*a)[0, (+/-)2/sqrt(3)].
They are expressed in terms of the supercell reciprocal vectors as follows:
(+/-)AA_inv (+/-)2/3*BB_inv
and
(+/-)4/3*BB_inv .
In other words, the K points are mapped onto (+/-)1/3*BB_inv.
So if your k-lattice division includes (0, 1/3) without shift
it will for sure hit a K-point.

Now, a solution to your problem:
- Choose the number of divisions along BB_inv
[the third line in the setting chosen, as the graphene plane is (y,z)]
divisible by 3;
- To get an isotrope k-sampling in the (y,z) plane,
choose the number of divisions along AA_inv [the third line]
larger than that along BB_inv - not 2 times larger than in your case,
just ~2/sqrt(3);
- add no shift.

Good luck.

Andrei Postnikov



> Dear all,
>
> I would be thankful if you help me with this question. Namely, I am
> studying a graphene TranSIESTA example located at
> http://dipc.ehu.es/frederiksen/tstutorial/index.php/Graphene. A
> tarball is provided there which contains two folders to calculate the
> transmission of graphene sheet. My question is regarding the
> calculations in the "small-supercell" folder. Particularly, after the
> scattering calculation is finished, when I use tbtrans to calculate
> the transmission, if I change the Monkhorst-Pack sampling in the .fdf
> file to the following:
>
> %block kgrid_Monkhorst_Pack
>1   00  0.0
>0   20  0  0.0
>0   0   10 0.0
> %endblock Kgrid_Monkhorst_Pack
>
> the transmission "does not come to zero at zero energy". But "it comes
> to zero when using the original Monkhorst-Pack grid", i.e. :
>
> %block kgrid_Monkhorst_Pack
>1   00  0.0
>0   20  0  0.5
>0   0   10 0.5
> %endblock Kgrid_Monkhorst_Pack
>
> As can be seen, I have only changed the values of displacements in the
> %block kgrid_Monkhorst_Pack and nothing else was changed in the .fdf
> file. Apparently, by changing the displacements, I have missed the K
> point of graphene but I can not figure out the mechanism of this. Any
> help to clarify what is going on is indeed welcome.
>
> Bests,
> Mohammad,
>

Re: [SIESTA-L] how to use AtomicCoordinatesOrigin

[apostnik]

Dear Chiu Fang Yu :

1. The block AtomicCoordinatesOrigin has 3 entries,
not 3*(the number of your test calculations).
In your example, I'd guess that Siesta reads the first three values
and goes along with it.
In order to try N values, you need to run N different calculations.

2. To be sure that Siesta digested the values you passed
the way you intended, check the output file (that is useful anyway;
nobody is so perfect to rule out possible misprints of keywords
in the input file), or - in what concerns the coordinates - the .XV file.

Best regards

Andrei Postnikov

> Dear:
>
> I want to use AtomicCoordinatesOrigin to check egg box effect.
> I typed the parameter like below:
> %block AtomicCoordinatesOrigin
> 0.000 0.000 0.000
> 0.000 0.000 0.007
> 0.000 0.000 0.014
> 0.000 0.000 0.020
> 0.000 0.000 0.027
> 0.000 0.000 0.034
> 0.000 0.000 0.041
> 0.000 0.000 0.050
> 0.000 0.000 0.055
> 0.000 0.000 0.062
> 0.000 0.000 0.068
> 0.000 0.000 0.075
> 0.000 0.000 0.082
> 0.000 0.000 0.089
> 0.000 0.000 0.096
> 0.000 0.000 0.103
> 0.000 0.000 0.109
> %endblock AtomicCoordinatesOrigin
>
> 0.007 is set by the reason that my system in z direction = 2.624, and when
> I set MeshCutoff = 200 Ry, the number of grid in z direction is 24, we set
> 16 steps in z direction, then we have 2.624/24/16=0.007
>
> I want to  get the force and energy from every step.
>
> but the siesta calculation just behave like I did not set the block
> AtomicCoordinatesOrigin
>
> I have set MD.UseSaveXV = F
>
> is there any other parameter that I should set or how do I submit the
> calculation job to get the force and energy vs every step to check egg
> box?
>
> Any suggestion is welcome.
> Thanks!
> --
> Best Regards
>
> 邱芳瑜  Chiu Fang Yu
>

Re: [SIESTA-L] chain relaxation

[apostnik]

The basis set a priori has nothing to do with symmetry.
You want to run a chain geometry; this means that the atoms
in the (x,y) plane are far enough from each other that
their basis functions do not overlap. Consequently,
there is no force nor stress in this direction;
whatever comes out is numerical noise.
Why don't you simply fix the lattice instead of
using variable cell.
However, fixing (correctly chosen) stress components
may effectively serve the same purpose.
It is too difficult for me to figure out which components
to fix...

Good luck

Andrei

andrei.postni...@univ-lorraine.fr


> I fix the stress to prevent it become the non-symmetry structure.
> since siesta is a software using localized basis sets not plane wave basis
> set,
> if I didn't stress the angle, it will become a non-symmetry structure.
>
>
>
> --
> Best Regards
>
> 邱芳瑜  Chiu Fang Yu
> 國立成功大學 材料科學與工程學系碩二
> MOBILE：0930287221(中華)
> GMAIL：joyce79...@gmail.com 
>
>
>
> 2014-12-17 14:56 GMT+08:00 joyce79928cc . :
>>
>> sorry, it is​
>> 'change the kpoint until the energy/atom (eV) change is less than 0.01eV
>> ?​'
>>
>> --
>> Best Regards
>>
>> 邱芳瑜  Chiu Fang Yu
>> 國立成功大學 材料科學與工程學系碩二
>> MOBILE：0930287221(中華)
>> GMAIL：joyce79...@gmail.com 
>>
>>
>>
>> 2014-12-17 14:55 GMT+08:00 joyce79928cc . :
>>>
>>> I was wondering how to decide whether it is converge?
>>> change the kpoint until the energy(eV) change is less than 0.01eV
>>> or...?
>>> Thanks!
>>>
>>> --
>>> Best Regards
>>>
>>> 邱芳瑜  Chiu Fang Yu
>>> 國立成功大學 材料科學與工程學系碩二
>>> MOBILE：0930287221(中華)
>>> GMAIL：joyce79...@gmail.com 
>>>
>>>
>>>
>>> 2014-12-17 14:51 GMT+08:00 joyce79928cc . :

 Thanks for reply. I will try to fix it!

 --
 Best Regards

 邱芳瑜  Chiu Fang Yu
 國立成功大學 材料科學與工程學系碩二
 MOBILE：0930287221(中華)
 GMAIL：joyce79...@gmail.com 



 2014-12-17 14:38 GMT+08:00 :
>
> Dear joyce79928:
>
> just some ideas.
> Your results may be noise due to insufficient basis and/or
> insufficient k-mesh ans/or insufficient MeshCutoff
> (probably all three).
> Does your DOS have any resemblance with what is expected for a 1-dim.
> chain?
> Does your stress definition in the GeometryConstrains block
> serve any purpose?
> There is a number of works with Siesta on Au chains around,
> maybe scroll through them first.
> Anyway, you can systematically proceed in the following way:
>
> Take just one atom along the chain step, not three
> (make the third lattice vector shorter :-).
> Fix it.
> For the k-mesh, you need just 1 point in the plane but more than 4
> along Z.
> Run a sequence of calculations just varying the 3d lattice vector
> (chain step) - single point, no relaxation.
> Trace the total energy curve as function of chain step.
> This will give you an approximation to "true Siesta" result.
> Study how it depends on basis, k-mesh and MeshCutoff.
> Find "good" values for all three.
> Make the usual "egg box" test.
> Then hopefully your forces will be OK so that
> you can add more atoms to your unit, release them, let them relax
> etc.
>
> Good luck
>
> Andrei
>
> andrei.postni...@univ-lorraine.fr
>
>
> > Dear:
> >
> > I want to relax the chain with Siesta.
> > but different atom position with same parameter will get different
> > structure.
> > for example, I set the distance between chain to 2.93 A and get a
> final
> > structure of 2.569 A between chain.
> > Then I change initial chain distance to 2.88 A and get a final
> structure
> > of
> > 2.758 A between chain.
> >
> > why this happened?
> >
> > Thanks for reading my mail and it is welcome to discuss with me.
> >
> > below's are my fdf input file.
> >
> > SystemName  bulk_au
> > SystemLabel bulk_au
> >
> > ==
> > ==
> > # SPECIES AND BASIS
> >
> > # Number of species
> > NumberOfSpecies 1
> > %block ChemicalSpeciesLabel
> >   1  79 Au
> > %endblock ChemicalSpeciesLabel
> >
> > PAO.BasisSizeSZP
> > PAO.EnergyShift  0.04 Ry
> >
> >
> > # K-points
> >
> > %block kgrid_Monkhorst_Pack
> > 4   0   0   0.0
> > 0   4   0   0.0
> > 0   0   4   0.5
> > %endblock kgrid_Monkhorst_Pack
> >
> >
> > # UNIT CELL AND ATOMIC POSITIONS
> >
> > # UNIT CELL
> > LatticeConstant   1.0 Ang
> > %block LatticeVectors
> > 8.81553253 0. 0.
> > 4.40776625 7.63447409 0.
> > 0. 0. 8.79434817
> > %endblock LatticeVectors
> >
> > # Atomic coordinates
> > NumberOfAtoms 3
> > AtomicCoordinatesFormat ScaledCartesian
> > %block AtomicCoordinatesAndAtomicSpecies
> > 4.407766232 4.241374236 0.

Re: [SIESTA-L] chain relaxation

[apostnik]

Dear Chiu Fang Yu:

The criterion of convergence varies from case to case.
It depends on which property you are after.
In the test I suggest, this is the minimum in the Etot(z-step) curve.
You want to get your chain step with a certain accuracy.
Each set of parameters will give you a different curve.
Enhance the parameters till the minimum of the curve
(also its curvature) gets stable.
The absolute value of energies in each curve is not important,
it may drift with parameters.

Best regards

Andrei

andrei.postni...@univ-lorraine.fr


> sorry, it is​
> 'change the kpoint until the energy/atom (eV) change is less than 0.01eV
> ?​'
>
> --
> Best Regards
>
> 邱芳瑜  Chiu Fang Yu
> 國立成功大學 材料科學與工程學系碩二
> MOBILE：0930287221(中華)
> GMAIL：joyce79...@gmail.com 
>
>
>
> 2014-12-17 14:55 GMT+08:00 joyce79928cc . :
>>
>> I was wondering how to decide whether it is converge?
>> change the kpoint until the energy(eV) change is less than 0.01eV or...?
>> Thanks!
>>
>> --
>> Best Regards
>>
>> 邱芳瑜  Chiu Fang Yu
>> 國立成功大學 材料科學與工程學系碩二
>> MOBILE：0930287221(中華)
>> GMAIL：joyce79...@gmail.com 
>>
>>
>>
>> 2014-12-17 14:51 GMT+08:00 joyce79928cc . :
>>>
>>> Thanks for reply. I will try to fix it!
>>>
>>> --
>>> Best Regards
>>>
>>> 邱芳瑜  Chiu Fang Yu
>>> 國立成功大學 材料科學與工程學系碩二
>>> MOBILE：0930287221(中華)
>>> GMAIL：joyce79...@gmail.com 
>>>
>>>
>>>
>>> 2014-12-17 14:38 GMT+08:00 :

 Dear joyce79928:

 just some ideas.
 Your results may be noise due to insufficient basis and/or
 insufficient k-mesh ans/or insufficient MeshCutoff
 (probably all three).
 Does your DOS have any resemblance with what is expected for a 1-dim.
 chain?
 Does your stress definition in the GeometryConstrains block
 serve any purpose?
 There is a number of works with Siesta on Au chains around,
 maybe scroll through them first.
 Anyway, you can systematically proceed in the following way:

 Take just one atom along the chain step, not three
 (make the third lattice vector shorter :-).
 Fix it.
 For the k-mesh, you need just 1 point in the plane but more than 4
 along Z.
 Run a sequence of calculations just varying the 3d lattice vector
 (chain step) - single point, no relaxation.
 Trace the total energy curve as function of chain step.
 This will give you an approximation to "true Siesta" result.
 Study how it depends on basis, k-mesh and MeshCutoff.
 Find "good" values for all three.
 Make the usual "egg box" test.
 Then hopefully your forces will be OK so that
 you can add more atoms to your unit, release them, let them relax etc.

 Good luck

 Andrei

 andrei.postni...@univ-lorraine.fr


 > Dear:
 >
 > I want to relax the chain with Siesta.
 > but different atom position with same parameter will get different
 > structure.
 > for example, I set the distance between chain to 2.93 A and get a
 final
 > structure of 2.569 A between chain.
 > Then I change initial chain distance to 2.88 A and get a final
 structure
 > of
 > 2.758 A between chain.
 >
 > why this happened?
 >
 > Thanks for reading my mail and it is welcome to discuss with me.
 >
 > below's are my fdf input file.
 >
 > SystemName  bulk_au
 > SystemLabel bulk_au
 >
 > ==
 > ==
 > # SPECIES AND BASIS
 >
 > # Number of species
 > NumberOfSpecies 1
 > %block ChemicalSpeciesLabel
 >   1  79 Au
 > %endblock ChemicalSpeciesLabel
 >
 > PAO.BasisSizeSZP
 > PAO.EnergyShift  0.04 Ry
 >
 >
 > # K-points
 >
 > %block kgrid_Monkhorst_Pack
 > 4   0   0   0.0
 > 0   4   0   0.0
 > 0   0   4   0.5
 > %endblock kgrid_Monkhorst_Pack
 >
 >
 > # UNIT CELL AND ATOMIC POSITIONS
 >
 > # UNIT CELL
 > LatticeConstant   1.0 Ang
 > %block LatticeVectors
 > 8.81553253 0. 0.
 > 4.40776625 7.63447409 0.
 > 0. 0. 8.79434817
 > %endblock LatticeVectors
 >
 > # Atomic coordinates
 > NumberOfAtoms 3
 > AtomicCoordinatesFormat ScaledCartesian
 > %block AtomicCoordinatesAndAtomicSpecies
 > 4.407766232 4.241374236 0. 1
 > 4.407766232 4.241374236 2.93144939 1
 > 4.407766232 4.241374236 5.86289878 1
 > %endblock AtomicCoordinatesAndAtomicSpecies
 >
 > ==
 > ==
 > # General variables
 >
 > ElectronicTemperature  300 K
 > MeshCutoff   200. Ry
 > xc.functional LDA   # Exchange-correlation
 functional
 > xc.authorsCA
 > SpinPolarized .false.
 > SolutionMethod Diagon
 >
 >
 > # SCF variables
 >
 > DM.MixSCF1   T
>

Re: [SIESTA-L] chain relaxation

[apostnik]

Dear joyce79928:

just some ideas.
Your results may be noise due to insufficient basis and/or
insufficient k-mesh ans/or insufficient MeshCutoff
(probably all three).
Does your DOS have any resemblance with what is expected for a 1-dim. chain?
Does your stress definition in the GeometryConstrains block
serve any purpose?
There is a number of works with Siesta on Au chains around,
maybe scroll through them first.
Anyway, you can systematically proceed in the following way:

Take just one atom along the chain step, not three
(make the third lattice vector shorter :-).
Fix it.
For the k-mesh, you need just 1 point in the plane but more than 4
along Z.
Run a sequence of calculations just varying the 3d lattice vector
(chain step) - single point, no relaxation.
Trace the total energy curve as function of chain step.
This will give you an approximation to "true Siesta" result.
Study how it depends on basis, k-mesh and MeshCutoff.
Find "good" values for all three.
Make the usual "egg box" test.
Then hopefully your forces will be OK so that
you can add more atoms to your unit, release them, let them relax etc.

Good luck

Andrei

andrei.postni...@univ-lorraine.fr


> Dear:
>
> I want to relax the chain with Siesta.
> but different atom position with same parameter will get different
> structure.
> for example, I set the distance between chain to 2.93 A and get a final
> structure of 2.569 A between chain.
> Then I change initial chain distance to 2.88 A and get a final structure
> of
> 2.758 A between chain.
>
> why this happened?
>
> Thanks for reading my mail and it is welcome to discuss with me.
>
> below's are my fdf input file.
>
> SystemName  bulk_au
> SystemLabel bulk_au
>
> ==
> ==
> # SPECIES AND BASIS
>
> # Number of species
> NumberOfSpecies 1
> %block ChemicalSpeciesLabel
>   1  79 Au
> %endblock ChemicalSpeciesLabel
>
> PAO.BasisSizeSZP
> PAO.EnergyShift  0.04 Ry
>
>
> # K-points
>
> %block kgrid_Monkhorst_Pack
> 4   0   0   0.0
> 0   4   0   0.0
> 0   0   4   0.5
> %endblock kgrid_Monkhorst_Pack
>
>
> # UNIT CELL AND ATOMIC POSITIONS
>
> # UNIT CELL
> LatticeConstant   1.0 Ang
> %block LatticeVectors
> 8.81553253 0. 0.
> 4.40776625 7.63447409 0.
> 0. 0. 8.79434817
> %endblock LatticeVectors
>
> # Atomic coordinates
> NumberOfAtoms 3
> AtomicCoordinatesFormat ScaledCartesian
> %block AtomicCoordinatesAndAtomicSpecies
> 4.407766232 4.241374236 0. 1
> 4.407766232 4.241374236 2.93144939 1
> 4.407766232 4.241374236 5.86289878 1
> %endblock AtomicCoordinatesAndAtomicSpecies
>
> ==
> ==
> # General variables
>
> ElectronicTemperature  300 K
> MeshCutoff   200. Ry
> xc.functional LDA   # Exchange-correlation functional
> xc.authorsCA
> SpinPolarized .false.
> SolutionMethod Diagon
>
>
> # SCF variables
>
> DM.MixSCF1   T
> MaxSCFIterations  300   # Maximum number of SCF iter
> DM.MixingWeight   0.03  # New DM amount for next SCF cycle
> DM.Tolerance  1.d-4 # Tolerance in maximum difference
> DM.UseSaveDM  true  # to use continuation files
> DM.NumberPulay 5
> #Diag.ParallelOverKyes
>
> # MD variables
>
> MD.FinalTimeStep 1
> MD.TypeOfRun CG
> MD.NumCGsteps 50
> MD.VariableCell  T
> UseSaveData  T
>
> # Output variables
>
> WriteMullikenPop1
> WriteBands  .false.
> SaveRho .false.
> SaveDeltaRho.false.
> SaveHS  .false.
> SaveElectrostaticPotential  True
> SaveTotalPotential  no
> WriteCoorXmol   .true.
> WriteMDXmol .true.
> WriteMDhistory  .false.
> WriteEigenvaluesyes
>
> %block GeometryConstraints
> stress 4 5 6
> position 1
> %endblock GeometryConstraints
>
>
>
> --
> Best Regards
>
> 邱芳瑜  Chiu Fang Yu
> 國立成功大學 材料科學與工程學系碩二
> MOBILE：0930287221(中華)
> GMAIL：joyce79...@gmail.com 
>

[SIESTA-L] tool to calculate Velocity AutoCorrelation Function after MD in Siesta

[apostnik]

To the attention of Siesta users interested in molecular dynamics:

I updated my homemade tool
to construct the velocity autocorrelation function from the files
created by a molecular dynamics run, i.e., .MD / .MD_CAR / .ANI ;
you can download from

http://www.home.uni-osnabrueck.de/apostnik/Software/VeloCorrF.zip

a small directory including Fortran sources and a simple Makefile,
to be edited to your compiling parameters. A documentation file is included
along with the distribution, which also covers some more general issues
of molecular dynamics with Siesta.

The result is just the density of modes
(velocity autocorrelation function Fourier-transformed and taken to square),
but, once the function is there, more fancy combinations can be programmed.

Enjoy, and feel free to share with me your encountered problems.

Good luck

Andrei

-- prof. Andrei Postnikov ---  andrei.postni...@univ-lorraine.fr
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
Inst. de Chimie, Physique et Matériaux, 1 Bd Arago, F-57078 Metz, France

Re: [SIESTA-L] Really Need Help on "NPT Ensemble"

[apostnik]

Dear Pedram,

I think it could be good if you carefully think about
what actually you want to achieve,
and why not get something to read about molecular dynamics.

You have a periodic system with 4 atoms per cell.
You say you relaxed your system in CG scheme; fine.
Now, molecular dynamics needs some kind of ensemble.
It all is about how the atoms wildly move in all different directions.
If you'd have a molecule, these 4 atoms would be all you have,
so running a MD simulation on them would be legitimate.
However, you impose periodic boundary conditions,
so the translated atoms in all cells are forced to move the same way.
This is not the way the Nature does it.
In other words, you want to explore the phase space
(of all particles in a crystal) yet constrain yourself
to an infinitesimal specially selected part of the phase space.
I'd not expect anything reasonable from such simulation
even if done without technical errors.
If you check the literature you'll see that MD simulations
are done either on molecules, or - in dense systems -
on large enough supercells.

So if you want to simulate how ZnO vibrates at some temperature,
you have to go to a much much larger supercell, no way around.

Good luck

Andrei Postnikov



> Hi,
>
> I really hope someone would kindly help me solve this problem.
>
> Enclosed, there's a .fdf file for bulk zno running for optimization under
> NPT ensemble.
> The geometry has been first relaxed up to 0.001 eV/Ang in CG scheme and
> then I started running it in MD.
>
> The system is not showing any convergence in step 1800 and the trends
> seems
> to be really odd.
>
> I would be grateful if you kindly take a look at the files to find the
> source of the problem.
>
> Best Regards,
> Pedram
>

Re: [SIESTA-L] Really Need Help on "NPT Ensemble"

[apostnik]

Dear Pedram,

I think it could be good if you carefully think about
what actually you want to achieve,
and why not get something to read about molecular dynamics.

You have a periodic system with 4 atoms per cell.
You say you relaxed your system in CG scheme; fine.
Now, molecular dynamics needs some kind of ensemble.
It all is about how the atoms wildly move in all different directions.
If you'd have a molecule, these 4 atoms would be all you have,
so running a MD simulation on them would be legitimate.
However, you impose periodic boundary conditions,
so the translated atoms in all cells are forced to move the same way.
This is not the way the Nature does it.
In other words, you want to explore the phase space
(of all particles in a crystal) yet constrain yourself
to an infinitesimal specially selected part of the phase space.
I'd not expect anything reasonable from such simulation
even if done without technical errors.
If you check the literature you'll see that MD simulations
are done either on molecules, or - in dense systems -
on large enough supercells.

So if you want to simulate how ZnO vibrates at some temperature,
you have to go to a much much larger supercell, no way around.

Good luck

Andrei Postnikov



> Hi,
>
> I really hope someone would kindly help me solve this problem.
>
> Enclosed, there's a .fdf file for bulk zno running for optimization under
> NPT ensemble.
> The geometry has been first relaxed up to 0.001 eV/Ang in CG scheme and
> then I started running it in MD.
>
> The system is not showing any convergence in step 1800 and the trends
> seems
> to be really odd.
>
> I would be grateful if you kindly take a look at the files to find the
> source of the problem.
>
> Best Regards,
> Pedram
>

RE: [SIESTA-L] reciprocal lattice vectors in SIESTA

[apostnik]

Dear Diana,
I don't think there is a rotation of nominal reciprocal vectors.
However, in case of doubt, you can generate the k-points
which are printed, and will give a clue. For example,
(1 1 1) divisions with 0.5 shift each
yield the half-endpoints of each reciprocal vector.

Best regards

Andrei Postnikov


> Thanks Abraham,
>
> My problem is more that I don't know the rotation angle. In VASP I often
> have the problem that the reciprocal lattice vectors I generate from the
> cross products are off by a rotation angle from the ones VASP generates. I
> will try to do it anyway.
>
> All the best,
>
> Diana
>
> Diana M. Otálvaro
> PhD Candidate
>
> Computational Material Science
> MESA+ Institute of Nanotechnology
> University of Twente.
> Enschede, Nederland
> 
> From: siesta-l-requ...@uam.es [siesta-l-requ...@uam.es] on behalf of
> Abraham Hmiel [abehm...@gmail.com]
> Sent: Friday, September 27, 2013 7:05 PM
> To: Siesta, Self-Consistent DFT LCAO program, http://www.uam.es/siesta
> Subject: Re: [SIESTA-L] reciprocal lattice vectors in SIESTA
>
> Hello Diana,
>
> As far as I'm aware, SIESTA doesn't print the reciprocal lattice vectors
> by default and no keyword exists to print them. However, it is relatively
> easy to create a short script in python or Matlab or such a program that
> will calculate them by inputting the unit cell vectors and then taking
> cross products (http://www.lcst-cn.org/Solid%20State%20Physics/Ch22.html)
> or, if you're ambitious, you can modify the code to print them in a new
> file or something.
>
> Regards,
>
> Abraham Hmiel
> Katherine Belz Groves Fellow in Nanoscience
> Xue Group, SUNY College of Nanoscale Science and Engineering
> http://abehmiel.net/about
>
>

Re: [SIESTA-L] HOMO/LUMO - plotting

[apostnik]

> Dear All,
>
> Thank you for help.
> I have  checked the values of all grids for my HOMO level and,
> unfortunatelly, IMAG and PHASE are not approximatelly zero.
> IMAG: [-0.515770, +0.470730]
> MOD: [0,0.771080]
> PHASE:[-0.732770,2.40900]
> REAL:[-0.5320400,0.573190]
> What does it mean?

Dear Karolina,
I'd say, then there must be either a bug,
or an error in the SIESTA input
(not a Gamma point really used),
or a bug/error in the way Denchar transformed your data.
To eliminate these one by one,

1. Check your output files for whether, indeed, it was Gamma
used in the calculation;
2. Look at bare eigenvectors (prior to where the wave functions
are reconstructed). Are they real?
...

Sorry for not being more specific.

> In my input file I have:
> %block WaveFuncKPoints
>0.00   0.00  0.   from 150 to 153
> %endblock WaveFuncKPoints
> So this is definietly the Gamma point.

And what do you have in your output?

> My second question is about this +/- values for real. What is the
> physical meaning of those values (signs) if I consider i.e. HOMO
> level?

There is no physical meaning, just an option to show two surfaces
IN THE SAME PLOT, provided by XCrySDen. The fact that they are expected
to be +/- the same value is simply because that's what the users want
most often (e.g., to figure out the shape of an antibonding orbital).
If you want to plot two different (arbitrary) levels, this is
possible after tweaking the data file (apply a constant shift to all
data points).

Best regards

Andrei Postnikov


>
> Best regards,
> Karolina Milowska
>
>
> 2013/9/17, apost...@uni-osnabrueck.de :
>> Dear Karolina,
>>
>> 1. In Gamma, the eigenvectors are (usually) real;
>> please have a careful look at the generated data;
>> I bet your IMAG and PHASE are roughly zero;
>> this is REAL what you need.
>>
>> 2. Taking MOD instead would mean to lose the information about the sign,
>> an important property of a wave function.
>>
>> 3. XCrySDen allows to plot +/- isolevels in the same plot
>> as surfaces of different color. Once you are in the
>> "Isosurface/Property-plane Controls"
>> (where you define the isovalue), press the button
>> "Render +/- isovalue"
>> You can define the properties of both isosurfaces by choosing
>> "Set COLOR prameters".
>>
>> Hopefully this will help
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>>
>>> Dear SIESTA Users,
>>>
>>> I would like to plot HOMO and LUMO level of my molecule.
>>> Therefore:
>>> 1) I read from EIG file numbers which correspond to those levels.
>>> 2) I run Denchar and got all functions (for Gamma point). I picked up
>>> those, which correspond to my HOMO and LUMO.
>>>
>>> But I do not know which one should I plot. I got 4 for each spin and
>>> wavefunction: IMAG or PHASE or REAL or MOD.
>>> Should I plot MOD (is this |phi|^2  ?) ?
>>>
>>> How to plot those levels (HOMO and LUMO) in Xcrysden in one picture
>>> (for both spins?) Is it possible? How to do that?
>>>
>>> Can anybody help me?
>>>
>>> Best regards,
>>> Karolina Milowska
>>>
>>
>>
>

Re: [SIESTA-L] HOMO/LUMO - plotting

[apostnik]

Dear Karolina,

1. In Gamma, the eigenvectors are (usually) real;
please have a careful look at the generated data;
I bet your IMAG and PHASE are roughly zero;
this is REAL what you need.

2. Taking MOD instead would mean to lose the information about the sign,
an important property of a wave function.

3. XCrySDen allows to plot +/- isolevels in the same plot
as surfaces of different color. Once you are in the
"Isosurface/Property-plane Controls"
(where you define the isovalue), press the button
"Render +/- isovalue"
You can define the properties of both isosurfaces by choosing
"Set COLOR prameters".

Hopefully this will help

Best regards

Andrei Postnikov



> Dear SIESTA Users,
>
> I would like to plot HOMO and LUMO level of my molecule.
> Therefore:
> 1) I read from EIG file numbers which correspond to those levels.
> 2) I run Denchar and got all functions (for Gamma point). I picked up
> those, which correspond to my HOMO and LUMO.
>
> But I do not know which one should I plot. I got 4 for each spin and
> wavefunction: IMAG or PHASE or REAL or MOD.
> Should I plot MOD (is this |phi|^2  ?) ?
>
> How to plot those levels (HOMO and LUMO) in Xcrysden in one picture
> (for both spins?) Is it possible? How to do that?
>
> Can anybody help me?
>
> Best regards,
> Karolina Milowska
>

Re: Fwd: [SIESTA-L] Electronic temperature in Anneal MD

[apostnik]

> dr. Roberto Pasianot, dr. Andrei Posnikov,
>
> Thanks for the replies.I have a question about the following:
>
> *"You may wish to start calculation with its value of say 300 - 600 K
> but then take care to reduce it to a small enough value that
> the total energy does not vary with it any more."*
>
> Does this mean in CG type of run and MD, once I have optimized my
> structure
> with low band-gap and 300K smearing, I should reduce the "electronic
> temperature" to go on calculating other properties, such as polarization?

Dear Pedram:
yes, you got the idea.
You can look at ElectronicTemperature as at yet another "cutoff",
to be used with care and checked against. For example:
you can converge your calculation with one k-point, but before
extracting the final result you may (should) wish to increase this number
and look how the "meaningful" results are affected. Similarly with
ElectronicTemperature: its elevated value helps the convergence but
spoils the results. It has any effect at all if you have a metal
or narrow-band semiconductor, i.e., when occupied and unoccupied
states are close enough, within your smearing value. "Other properties"
are affected inasmuch as they depend on the integration over the
Brillouin zone. You ask about polarization; this normally means that
you have a substantial band gap, probably much broader than the
ElectronicTemperature you will use.
Summarizing: normally you calculate the polarization
(or, other "final" properties) on a well converged case,
when the convergence is no more an issue and you can (should)
set the ElectronicTemperature to very small value without any problem.

Best regards

Andrei Postnikov



> Regards,
> Pedram
>
>
> On Thu, Sep 5, 2013 at 7:48 PM, Roberto Pasianot
> wrote:
>
>> Hi Pedram,
>>
>> MD temperature is "reaĺ", thermodynamic temperature.
>> Electronic temperature is just a smearing parameter for the
>> Fermi surface, that easies the computation of integrals in
>> reciprocal space involving occupied states. You could use
>> the MP scheme instead of the FD one, thus loosing connection
>> with any statistical distribution. Or even choose zero, with
>> possible horrible effects on SCF convergence.
>>
>> Bye, bye,
>>
>> Roberto
>>
>>
>>
>> On 09/05/2013 11:41 AM, pedram heidari wrote:
>>
>>> Hi,
>>>
>>> I have problem with temperatures in MD Anneal in Si.fdf that is present
>>> in
>>> "tests" directory in siesta code. (I have enclosed the .fdf)
>>>
>>> There, I can see the electronic temperature is 300K and two other
>>> parameters as:
>>>
>>> MD.InitialTemperature 100 K
>>> MD.TargetTemperature 600 K
>>>
>>> To my understanding, the electronic temperature comes from the
>>> fermi-dirac
>>> distribution and deals with the valence "electrons" but MD temperatures
>>> come
>>> from Maxwell-Boltzmann distribution for "atoms" and should consider the
>>> role of
>>> the nuclei as well, if I am not mistaken.
>>>
>>> What is exactly the point of mentioning the temperature in two
>>> different
>>> distributions in one MD simulation?
>>>
>>> Thanks,
>>> Pedram
>>>
>>>
>>
>

Re: [SIESTA-L] Electronic temperature in Anneal MD

[apostnik]

Hi Pedram,
the two "temperatures" are totally decoupled and stay in no relation
one to the other.

The ElectronicTemperature a priori has nothing to do
with molecular dynamics; its effect is on electronic structure
iterations, irrespectively of whether the calculation is static
or dynamic. As explained in the manual, it controls the width
of the Fermi step function. Think of it as merely an artificial
smearing parameter, to smoothen the electronic structure iterations.
It is particularly useful if you have a metal or narrow band gap,
where there are many levels close to HOMO - LUMO.
Note that the exact value of total energy depends on this parameter
quite a lot.
You may wish to start calculation with its value of say 300 - 600 K
but then take care to reduce it to a small enough value that
the total energy does not vary with it any more.

The temperature definition(s) in the MD part controls the thermostat
and hence the "classical" behaviour of atoms, without directly
affecting the electrons.

Best regards

Andrei Postnikov

> Hi,
>
> I have problem with temperatures in MD Anneal in Si.fdf that is present in
> "tests" directory in siesta code. (I have enclosed the .fdf)
>
> There, I can see the electronic temperature is 300K and two other
> parameters as:
>
> MD.InitialTemperature 100 K
> MD.TargetTemperature 600 K
>
> To my understanding, the electronic temperature comes from the fermi-dirac
> distribution and deals with the valence "electrons" but MD temperatures
> come from Maxwell-Boltzmann distribution for "atoms" and should consider
> the role of the nuclei as well, if I am not mistaken.
>
> What is exactly the point of mentioning the temperature in two different
> distributions in one MD simulation?
>
> Thanks,
> Pedram
>

Re: [SIESTA-L] One more question. I do apologize.

[apostnik]

Dear Andrei:

The calculation of band structure is just an extra calculation
over selected k points. Technically it is independent on whether
your previous (self-consistent) calculation was done just with Gamma,
or with more k-points, be it a single cell or a supercell.

As you pass to a supercell, the BZ is reduced; you have more (folded) bands
which disperse less, because the k-path becomes shorter.

Take care of how you define the path; there is a potential danger
to mess up absolute / relative units, and relative with respect
to single cell lattice constant, or the supercell one.

Best regards

Andrei Postnkov


> Let's say i have computed supercell with Gamma point only. Later i decided
> to look for a bandstructure and specify path in bandLines block G-X-M-L
> 
>
> I receive flat band picture, although it should not be the case. My
> question is, will Siesta sample at this points defined along different
> paths or is it completely wrong to start out with Gamma point
> calculations,
> and then specify BandLines block ?
>
> With Best Regards, Andrei Buin.
>

Re: [SIESTA-L] One more question. I do apologize.

[apostnik]

Dear Andrei:

The calculation of band structure is just an extra calculation
over selected k points. Technically it is independent on whether
your previous (self-consistent) calculation was done just with Gamma,
or with more k-points, be it a single cell or a supercell.

As you pass to a supercell, the BZ is reduced; you have more (folded) bands
which disperse less, because the k-path becomes shorter.

Take care of how you define the path; there is a potential danger
to mess up absolute / relative units, and relative with respect
to single cell lattice constant, or the supercell one.

Best regards

Andrei Postnkov


> Let's say i have computed supercell with Gamma point only. Later i decided
> to look for a bandstructure and specify path in bandLines block G-X-M-L
> 
>
> I receive flat band picture, although it should not be the case. My
> question is, will Siesta sample at this points defined along different
> paths or is it completely wrong to start out with Gamma point
> calculations,
> and then specify BandLines block ?
>
> With Best Regards, Andrei Buin.
>

Re: [SIESTA-L] Bizarre behaviour

[apostnik]

OK Andrei.
Now the chemical formula makes sense.
(Is the methyl ammonium rotating in its site?
Disordered in this sense? - Must be funny... - A nice task for Siesta)

1. The 2x2x2 supercell would nicely map X, M and R points
of the primitive cell onto its new Gamma. Hopefully it's easy
to trace "by hand", checking contributions from different atoms
to the wavefunction at these points in case of difficulty.
I don't understand why do you need a 3x3x3 supercell if you are
interested just in folding. Must be a long calculation...

2. Be careful with the supercell block; my suggestion is
to use it just for generating the supercell, then put the generated
coordinates into the input (or start from XV) and deactivate
the supercell block. Otherwise it is easy to make an error.

Good luck

Andrei Postnikov


> Thank you Andrei.
>
> 1. I do apologize for the typo. System is pervoskite (CH3NH3)SnI3. It has
> 12 atoms in cubic cell.
> 2. Going to supercell should not increase the bandgap. I agree on this.
> 3. I did tetragonal cell of the same system of 48 atoms , and it looks ok
> if i  use superecell but it has bandgap  at G-point.
> 4. I looked structures before, basically 12 atoms at 8x8x8 k-mesh gives
> bandgap of 0.95 eV, whereas with %blocksupercell 3 3 3 %endblcoksupercell
> with only G point  gives 1.6 eV.
> 5. And you are right, i also started to suspect that zone folding is the
> problem(i mean, my whole point of doing supercell calculations is to map
> R-point into G-point).
> 6. Visually it looks ok. It contains methylammonium. The only thing it
> might have is the dipole moment.
>
> Witll try to use 2x2x2,4x4x4 (It should correctly map R point onto G
> point,
> correct me if i'm wrong )  supercell to see whether problem disappears.
>
>
> With Best Regards, Andrei Buin.
>
>
>
> On Wed, Aug 21, 2013 at 6:29 AM,  wrote:
>
>> Dear Andrei:
>> I don't see how your system can be a perovskite,
>> nor how CH3N3SnI3 sums up to 12 atoms. I'll appreciate if you
>> give a hint, just for my education.
>> However, this is beyond the point.
>> In principle, a mere passing to a supercell
>> should not result in an increase of the band gap.
>> The mapping of k-points may be not exact, i.e.,
>> Gamma of the 3*3*3 supercell maps onto 0, 1/3, 2/3 [*(2*pi/a)]
>> divisions along axes of the reciprocal cell,
>> none of which points includes R=(1/2 1/2 1/2). You can check
>> what the band gap (in the primitive cell) is at (1/3 1/3 1/3),
>> or at (3/8 3/8 3/8), the closest of your k-points at the 8x8x8 mesh.
>> [However, you say that you have an automatic k sampling, then
>> the chances are that your mesh is shifted... and does not contain
>> the R point..?]
>> Apart from these tests, I think that some simple error in structure
>> (as simple as a wrong lattice constant in a supercell, sorry)
>> is a very good candidate for this kind of problem.
>> To check it, it might be a good idea to have a view at both structures
>> (the primitive one and the supercell) as stored in corresponding XV
>> files,
>> and to compare their corresponding DOS.
>>
>> Again apart from all this, if your system contains methyl azide,
>> I think it could be tricky to initialize the different charge states
>> at different N atoms correctly... If you just let it go as it wish
>> the calculation may converge to something wrong... Even to differently
>> wrong in two different supercells.
>>
>> Sorry for too many suggestions.
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>> > Dear Siesta users,
>> >
>> > I'm having very bizarre problem.
>> > I'm simulating the primitive cubic unit cell(12 atoms) of the
>> CH3N3SnI3
>> > perovksite  with GGA and PBE and gives me band gap of Eg=0.9 eV at R
>> point
>> > with 8x8x8 Moshkrof automatic k sampling(i know it's wrong in DFT, but
>> > still) and then i use supercell approach built  in  with 3x3x3 unit
>> > cell(432 atoms) with Gamma point only and i get a bangap of 1.6 eV.
>> I've
>> > checked against plane cutoff convergence seems ok.  I thought maybe
>> > dispersion  is too big so one should sample it more finely, but it
>> turns
>> > out  meff=0.13me, which should be fine. Any ideas ?
>> >
>> >
>> > With Best Regards, Andrei Buin.
>> >
>>
>>
>

Re: [SIESTA-L] Bizarre behaviour

[apostnik]

Dear Andrei:
I don't see how your system can be a perovskite,
nor how CH3N3SnI3 sums up to 12 atoms. I'll appreciate if you
give a hint, just for my education.
However, this is beyond the point.
In principle, a mere passing to a supercell
should not result in an increase of the band gap.
The mapping of k-points may be not exact, i.e.,
Gamma of the 3*3*3 supercell maps onto 0, 1/3, 2/3 [*(2*pi/a)]
divisions along axes of the reciprocal cell,
none of which points includes R=(1/2 1/2 1/2). You can check
what the band gap (in the primitive cell) is at (1/3 1/3 1/3),
or at (3/8 3/8 3/8), the closest of your k-points at the 8x8x8 mesh.
[However, you say that you have an automatic k sampling, then
the chances are that your mesh is shifted... and does not contain
the R point..?]
Apart from these tests, I think that some simple error in structure
(as simple as a wrong lattice constant in a supercell, sorry)
is a very good candidate for this kind of problem.
To check it, it might be a good idea to have a view at both structures
(the primitive one and the supercell) as stored in corresponding XV files,
and to compare their corresponding DOS.

Again apart from all this, if your system contains methyl azide,
I think it could be tricky to initialize the different charge states
at different N atoms correctly... If you just let it go as it wish
the calculation may converge to something wrong... Even to differently
wrong in two different supercells.

Sorry for too many suggestions.

Best regards

Andrei Postnikov


> Dear Siesta users,
>
> I'm having very bizarre problem.
> I'm simulating the primitive cubic unit cell(12 atoms) of the CH3N3SnI3
> perovksite  with GGA and PBE and gives me band gap of Eg=0.9 eV at R point
> with 8x8x8 Moshkrof automatic k sampling(i know it's wrong in DFT, but
> still) and then i use supercell approach built  in  with 3x3x3 unit
> cell(432 atoms) with Gamma point only and i get a bangap of 1.6 eV. I've
> checked against plane cutoff convergence seems ok.  I thought maybe
> dispersion  is too big so one should sample it more finely, but it turns
> out  meff=0.13me, which should be fine. Any ideas ?
>
>
> With Best Regards, Andrei Buin.
>

Re: [SIESTA-L] PDOS: Is there a way to specify radius of projection?

[apostnik]

Hi Benedikt,
Alberto is right, SIESTA does not need projections,
but YOU might need them; exactly for the sake of comparing
the charges with those from "muffin-tin"-type codes
I wrote some time ago a primitive tool
http://www.home.uni-osnabrueck.de/apostnik/Software/grdint.f
which "integrates" grid properties, e.g. charge (spin) densities,
over given (atom-centered) spheres.
Some remarks:
1. It integrates functions defined on the grid, which are not
(l,m) resolved. That means, you can produce spin-up and spin-down
charges, but not partial s,p,d-charges. Of course you can first
generate LDOS within some energy interval, if you find it useful,
and then integrate it over a sphere.
2. The "integration" is in fact merely a counting of grid points
which either fall within, or not, of a given sphere. So the result
is not very accurate and prone to "noise". But it is usually OK
for the sake of comparison, and anyway becomes "better"
as the mesh density is increased.
3. The tool as not fast as you may expect it to be, for the task
it performs, because the algorithm used is very straightforward;
please feel free to improve.

Best regards

Andrei Postnikov



> Hi Benedikt,
>
> SIESTA does not need projections, as the basis orbitals are localized on
> the atoms. You get naturally the "chemical" information one is used to.
> Plane-wave codes such as vasp do need projections to get some kind of
> "local" information from the delocalized basis.
>
>   Alberto
>
>
> On Mon, Jun 17, 2013 at 4:28 PM, Benedikt Ziebarth <
> benedikt.zieba...@kit.edu> wrote:
>
>> Hello,
>>
>> I have a question about PDOS calculations with siesta. Is there a way to
>> specify the radius around the atoms in which the projection is carried
>> out?
>> This option exists in different other dft codes like vasp (
>> http://cms.mpi.univie.ac.at/**vasp/vasp/RWIGS.html<http://cms.mpi.univie.ac.at/vasp/vasp/RWIGS.html>
>> ).
>> Any help would be very welcome.
>> Thanks
>>
>> Benedikt Ziebarth
>>
>

Re: [SIESTA-L] Problems with equilibrium lattice constant

[apostnik]

Dear Andreas,
I'd suggest you first to check whether the band structure,
and not just the band gap, comes out right.
Try to have a look at least at PDOS.
The point is, NiMnSb has some magnetic structure,
which would hardly emerge by itself by default.
You define spin polarization = true in your input file
but do not elaborate further. In this case (see the Siesta manual),
all atoms (including Sb) are initialized with maximal magnetic moment,
all parallel. The chances are, such initial configuration
will converge, but to something completely unreasonable.
You can get an idea of the resulting magnetic structure
by switching on mulliken occupations,
Probably it converges to different magnetic structures at
different lattice constant, hence the "ugly jumps".
Another possible reason for discontinuities in the Etot(a) curve
may be a too low k-grid cutoff, so that the actual (low)
k-mesh number of divisions changes at some values of a.

Best regards

Andrei Postnikov

>
> I'm a new user to siesta and want to determine the equilibrium
> lattice constant for the half-metallic heusler alloy NiMnSb. To do this I
> downloaded the GGA pseudopotentials for Ni and Mn and created my own for
> Sb with the atom program. My problem is that either I have no band gap at
> the fermi Energy in the minority spin domain , or the plot total Energy
> vs. lattice constant has ugly jumps.  I could not yet achieve a nice
> parabolic curve for the lattice constant with a band gap at the minimum
> Value. I'll attach my input file. If any other files are important to
> you just let me know. I'd be glad if anyone can help me out.
>
>  
>
> greetings, Andreas Prinz-Zwick

Re: [SIESTA-L] Convergence for charged system...

[apostnik]

Dear Luigi,

Si5H12 has an even number of electrons, therefore
when adding or removing one, in principle, you ought to use
spin polarized = true.
(This will become less and less crucial as you increase
the system so it becomes more "solid-like" and/or metallic.
But at the beginning this extra electron is likely to be
carried by a isolated level in the gap, I guess).
A bad convergence is probably due to swapping between
competing states, where this extra electron is differently
allocated in space. This might be a manifestation of
a genuine physical problem in a sense that,
given the symmetry of your system, the level in question
may want to be degenerate (jus a guess...) whereas you
populate it with a single electron. If this is indeed the case,
a small distortion may help to fix the issue... However
"in real world" you'd need to allow spin-orbit interaction
to get it right.

By setting spin polarization off you cheat the system and tell it
to search only among doubly degenerate states. The result may
converge better, but place the energy level at a wrong place.
It could be useful to have a look at where the levels are
(where your extra charge is sitting) and try to make sense
out of it.

If not caring about all above, but just technically
aiming at convergence, I'd suggest you to
i) reduce the mixing factor considerably,
ii) first set ElectronicTemperature higher (say to 600 K)
to facilitate convergence, and then gradually decrease it.

Moreover the MeshCutoff of 150 Ry might be not high enough
to get a good physical result, and might also affect the convergence.

Good luck

Andrei Postnikov


> Hello everybody,
> I am trying to perform a calculation of a charged system composed of a
> molecule of Si5H12
> with an extra electron (and in the future also one less electron) so I
> set:
>
> NetCharge -1.00
> SpinPolarized .true.
> FixSpin   .true.
> TotalSpin 1.0
>
> But the calculation even for this small system does not converge.
> If I put SpinPolarization .false. the calculation converges rapidly but I
> am not sure that it's correct.
> recently I added also:
>
> SimulateDoping   .true.
>
> It helps only slightly but convergence is still far from beying achieved.
> These are the other parameters of my input file:
> --
> SystemName   SiH
> SystemLabel   SiH
> #(1) SPECIES AND ATOMS
> NumberOfAtoms  17
> NumberOfSpecies 2
> %block ChemicalSpeciesLabel
>  1  14  Si # Species index, atomic number, species label
>  2   1   H # Species index, atomic number, species label
> %endblock ChemicalSpeciesLabel
> #(2) FUNCTIONAL
> xc.functional  GGA
> xc.authors PBE
> SpinPolarized .true.
> FixSpin   .true.
> TotalSpin 1.0
> MeshCutoff150 Ry
> NetCharge -1.00
> SimulateDoping   .true.
> #(3) BASIS
> PAO.SplitTailNorm   .true.
> #(4) SCF OPTIONS
> MaxSCFIterations  4000
> DM.MixingWeight0.01
> DM.Tolerance   1.d-4
> DM.NumberPulay 3
> SolutionMethod Diagon
> ElectronicTemperature   300 K
> #(5) MD OPTIONS
> MD.TypeOfRun   CG
> MD.NumCGsteps  1000
> #(6) OUTPUT OPTIONS
> #(7) CRYSTAL STRUCTURE
> LatticeConstant 30.0 Ang
> %block LatticeVectors
> 1.000 0.000 0.000
> 0.000 1.000 0.000
> 0.000 0.000 1.000
> %endblock LatticeVectors
> #(8) ATOMIC COORDINATES
> AtomicCoordinatesFormat  NotScaledCartesianAng
> AtomicCoordinatesAndAtomicSpecies < Si-H.siesta
> --
>
> I hope someone has a hint.
> Thanck you very much in advance.
>
> L.
>

Re: [SIESTA-L] Convergence for charged system...

[apostnik]

Dear Luigi,

Si5H12 has an even number of electrons, therefore
when adding or removing one, in principle, you ought to use
spin polarized = true.
(This will become less and less crucial as you increase
the system so it becomes more "solid-like" and/or metallic.
But at the beginning this extra electron is likely to be
carried by a isolated level in the gap, I guess).
A bad convergence is probably due to swapping between
competing states, where this extra electron is differently
allocated in space. This might be a manifestation of
a genuine physical problem in a sense that,
given the symmetry of your system, the level in question
may want to be degenerate (jus a guess...) whereas you
populate it with a single electron. If this is indeed the case,
a small distortion may help to fix the issue... However
"in real world" you'd need to allow spin-orbit interaction
to get it right.

By setting spin polarization off you cheat the system and tell it
to search only among doubly degenerate states. The result may
converge better, but place the energy level at a wrong place.
It could be useful to have a look at where the levels are
(where your extra charge is sitting) and try to make sense
out of it.

If not caring about all above, but just technically
aiming at convergence, I'd suggest you to
i) reduce the mixing factor considerably,
ii) first set ElectronicTemperature higher (say to 600 K)
to facilitate convergence, and then gradually decrease it.

Moreover the MeshCutoff of 150 Ry might be not high enough
to get a good physical result, and might also affect the convergence.

Good luck

Andrei Postnikov


> Hello everybody,
> I am trying to perform a calculation of a charged system composed of a
> molecule of Si5H12
> with an extra electron (and in the future also one less electron) so I
> set:
>
> NetCharge -1.00
> SpinPolarized .true.
> FixSpin   .true.
> TotalSpin 1.0
>
> But the calculation even for this small system does not converge.
> If I put SpinPolarization .false. the calculation converges rapidly but I
> am not sure that it's correct.
> recently I added also:
>
> SimulateDoping   .true.
>
> It helps only slightly but convergence is still far from beying achieved.
> These are the other parameters of my input file:
> --
> SystemName   SiH
> SystemLabel   SiH
> #(1) SPECIES AND ATOMS
> NumberOfAtoms  17
> NumberOfSpecies 2
> %block ChemicalSpeciesLabel
>  1  14  Si # Species index, atomic number, species label
>  2   1   H # Species index, atomic number, species label
> %endblock ChemicalSpeciesLabel
> #(2) FUNCTIONAL
> xc.functional  GGA
> xc.authors PBE
> SpinPolarized .true.
> FixSpin   .true.
> TotalSpin 1.0
> MeshCutoff150 Ry
> NetCharge -1.00
> SimulateDoping   .true.
> #(3) BASIS
> PAO.SplitTailNorm   .true.
> #(4) SCF OPTIONS
> MaxSCFIterations  4000
> DM.MixingWeight0.01
> DM.Tolerance   1.d-4
> DM.NumberPulay 3
> SolutionMethod Diagon
> ElectronicTemperature   300 K
> #(5) MD OPTIONS
> MD.TypeOfRun   CG
> MD.NumCGsteps  1000
> #(6) OUTPUT OPTIONS
> #(7) CRYSTAL STRUCTURE
> LatticeConstant 30.0 Ang
> %block LatticeVectors
> 1.000 0.000 0.000
> 0.000 1.000 0.000
> 0.000 0.000 1.000
> %endblock LatticeVectors
> #(8) ATOMIC COORDINATES
> AtomicCoordinatesFormat  NotScaledCartesianAng
> AtomicCoordinatesAndAtomicSpecies < Si-H.siesta
> --
>
> I hope someone has a hint.
> Thanck you very much in advance.
>
> L.
>

Re: [SIESTA-L] phonon graphene

[apostnik]

> Dear siesta user
>
> I got graphene phonon dispersion with vibra. without
> kgrid_Monkhorst_Pack(%block), means with just single k-point (gamma)
> and   superr cell 2*2*0  got to a accurate graphene dispersion.
>
> I think that it's not true physically but result is correct.
> my question is :
>  Could a single k-point along the  gamma direction be enough to
> describe phonon dispersion?
>
> Best
> Drogar
>

Dear Drogar,
it seems that you mix up two issues:

Monkhorst_Pack, single k-point or not define the quality
(accuracy) of your electronic structure calculation. I'd say,
single k-point for graphene is not good enough, but you'll certainly
get some levels. However, you'll have no dispersion of electron bands,
no Dirac point etc.

The phonon dispersion comes from how frequency varies with wavevector
of phonon, say q. It has nothing to do with Monkhorst-Pack;
technically in Siesta it is given by supercell size. In your case
supercell has more than one cell in the (X,Y) plane, so you ought
to get SOME phonon dispersion. The question of how good it is can be
answered by testing the results vs. increased supercell size
(I bet there is already enough publications dealing with this issue).
A priori 5*5*1 supercell is not obviously too small,
but it depends on which effects you are after.

Your "but result is correct" may mean that you've got luck and
the force constants are not too sensitive to your insufficient
k-mesh integration. Usually they are. Look more carefully;
maybe you results are not that correct after all.

Best regards

Andrei Postnikov

Re: [SIESTA-L] phonon graphene

[apostnik]

> Dear siesta user
>
> I got graphene phonon dispersion with vibra. without
> kgrid_Monkhorst_Pack(%block), means with just single k-point (gamma)
> and   superr cell 2*2*0  got to a accurate graphene dispersion.
>
> I think that it's not true physically but result is correct.
> my question is :
>  Could a single k-point along the  gamma direction be enough to
> describe phonon dispersion?
>
> Best
> Drogar
>

Dear Drogar,
it seems that you mix up two issues:

Monkhorst_Pack, single k-point or not define the quality
(accuracy) of your electronic structure calculation. I'd say,
single k-point for graphene is not good enough, but you'll certainly
get some levels. However, you'll have no dispersion of electron bands,
no Dirac point etc.

The phonon dispersion comes from how frequency varies with wavevector
of phonon, say q. It has nothing to do with Monkhorst-Pack;
technically in Siesta it is given by supercell size. In your case
supercell has more than one cell in the (X,Y) plane, so you ought
to get SOME phonon dispersion. The question of how good it is can be
answered by testing the results vs. increased supercell size
(I bet there is already enough publications dealing with this issue).
A priori 5*5*1 supercell is not obviously too small,
but it depends on which effects you are after.

Your "but result is correct" may mean that you've got luck and
the force constants are not too sensitive to your insufficient
k-mesh integration. Usually they are. Look more carefully;
maybe you results are not that correct after all.

Best regards

Andrei Postnikov

Re: [SIESTA-L] plot charge density and wave function

[apostnik]

Dear Somayeh Fotohi:

In XCrySDen, you select in menu
File -> Open Structure -> Gaussian98 Cube File

and read in your .cube

For the rest, check the XCrySDen documentation.

Good luck

Andrei Postnikov

>
> Dear siesta users,
>
>
> I m
> trying to plot charge density and wave function in  3D.
> For
> this I do:
> 1- Run
> siest.fdf for obtaining  the files required by  denchar program
> such as (. DM. DIM. PLD.
> WSF)  .ion of each species . 2- Run denchar.fdf program ( mode of
> representation ( 3D).)by  denchar denchar.fdf I obtain many file such as
> .PHASE.cube, REAL.cube, IMAG.cube, .DRHO.cube, .RHO.cube, …. but I dont
> know how to visualize by xcrysden correctly. Please help me.  I look
> forward to your kindly reply and suggestion. Best Regards, somayeh fotoohi

Re: [SIESTA-L] plot charge density and wave function

[apostnik]

Dear Somayeh Fotohi:

In XCrySDen, you select in menu
File -> Open Structure -> Gaussian98 Cube File

and read in your .cube

For the rest, check the XCrySDen documentation.

Good luck

Andrei Postnikov

>
> Dear siesta users,
>
>
> I m
> trying to plot charge density and wave function in  3D.
> For
> this I do:
> 1- Run
> siest.fdf for obtaining  the files required by  denchar program
> such as (. DM. DIM. PLD.
> WSF)  .ion of each species . 2- Run denchar.fdf program ( mode of
> representation ( 3D).)by  denchar denchar.fdf I obtain many file such as
> .PHASE.cube, REAL.cube, IMAG.cube, .DRHO.cube, .RHO.cube, …. but I dont
> know how to visualize by xcrysden correctly. Please help me.  I look
> forward to your kindly reply and suggestion. Best Regards, somayeh fotoohi

Re: [SIESTA-L] LUMO and HOMO

[apostnik]

Dear Somayeh Fotohi:
your question is not clear enough.

> How can I got the figures of LUMO and HOMO ?

LUMO and HOMO are molecular orbitals. They can be identified
knowing the number of valence electrons (possibly taking into account
different occupations in spin-up and spin-down channels).

Their energies can be found in the .EIG file
where the energies of all orbitals appear in increased order,
first spin 1 then spin 2.

For plotting, you either pick up wavefunctions individually
by their numbers, or
you plot LDOS in a given energy interval (which may of course include
just one orbital, if you wish).

> I need to  plot in 2D LDOS of my system integrated between two energy
> levels(HOMO and LUMO)

If integrated BETWEEN HOMO and LUMO you'll normally get an empty
data file, because there are no states by definition.

Best regards

Andrei Postnikov


> Dear Siesta users,
> I need to  plot in 2D LDOS of my system integrated between two energy
> levels(HOMO and LUMO)
> How can I got the figures of LUMO and HOMO ?
> somayeh fotoohi

Re: [SIESTA-L] Siesta utilities

[apostnik]

> Dear Siesta users,
>
> I want to know how to convert eig2bxsf and rho2xsf?. I have executable
> files of both, please tell me the next step.
> Thank u.
>

Dear Manjeet,
if you have executable files of both, then hopefully you've got them
as part of the whole distribution, which also includes the README file.
In a nutshell: you start the script in the command line
(type its name followed by "Enter") and then type in answers
to the queries the script asks you.
Typically, you do it in a directory where your files are ( .XV, .RHO, .EIG
and others).
Especially in case of rho2xsf, make sure that the executable is compiled
on the same architecture as your SIESTA code, because otherwise there might
be incompatibility in the binary format for .RHO; this is the most often
source of trouble.

Good luck

Andrei Postnikov

Re: [SIESTA-L] run by siesta

[apostnik]

> dear siesta users
> i want to run Ar gas, can i do this by siesta?
> if i can, are there any software which i could find coordinate of atoms in
> Ar gas?
> thanks

Dear roz_gh,
you question might have amused some colleagues, including myself,
but there is a substance in it.

"Yes you can", but you should consider that Siesta manipulates
"short-tail" objects; so "most of the time" in a gas there will be
no overlap and no interaction between your atoms. But, maybe that's
just OK for a gas - the atoms will develop forces only on collisions.
Comparing to a planewave method, Siesta might be
quite efficient. But how really good?

I doubt you'll find  software which gives you positions in Ar gas,
but a priori, just to start molecular dynamics, you'd need to generate
coordinates of atoms in a box of chosen size,
taken at a concentration of your choice,
using some random numbers generator. (Then you'll run MD for a while
anyway, I guess). This might be faster to program yourself
than search for a code which does it.

However, the real question is, what do you want achieve
by "running a gas"? To get its electronic structure and density
of states? - Hm. - To check the Mendeleyev-Klapeyron equation? - Good luck.
In fact, this may be didactically interesting. (Did anyone did
something similar before?)
Anyway, I suggest that you have some clear idea of what you want to do
before you start calculations.

For the future, please:
- give a meaningful subject line, to facilitate the search over archives,
- sign your messages

Best regards

Andrei Postnikov

Re: [SIESTA-L] MD-trajectory

[apostnik]

> Hi,
>
>   Is there any software (viewer) is available to view the MD trajectory
> output of siesta.
>
> regards
> hakkim
>

Dear Hakkim
you may check Sies2xsf from
http://www.home.uni-osnabrueck.de/apostnik/download.html
which contains md2axsf, to generate an input file for XCrySDen
(for an animation).

Best regards

Andrei Postnikov

Re: [SIESTA-L] Need explanation of force-constant matrix or force-constant in cartesian coordinate

[apostnik]

Dear Bedamani,

FC is organized as follows (following the first title line):
- the external loop over atoms which undergo finite displacements, i.e.,
  from MD.FCfirst to MD.FClast (or all atoms in the cell by default);
- for each atom, 6 displacements are applied in the order: -X,X,-Y,Y,-Z,Z ;
- for each displacement, the forces induced on ALL atoms in the cell
  are written, line by line, divided by magnitude of displecement,
  three Cartesian coordinates in a line.
  What is actually written is -(force)/(displacement) in units of eV/(Ang^2)

Best regards

Andrei Postnikov


> Dear siesta users,
>
>   We are trying to use the force-constant output from SIESTA for a
> 2-atom
>  (2-dimensions) crystal (graphene) as input into our program
>  which requires a force-constant matrix in Cartesian
>  coordinates.  The required full matrix would be 6x6 ([basicaly 3n*3n, for
> a n atom system],
>  with elements at the head Atom1x, Atom1y, Atom1z,Atom2x, Atom2y, Atom2z.
>  For SIESTA output from graphene we are getting 72
>  total elements which do not seem to be arranged in this order.
>  Can you tell us what the elements in the .FC file are, and how we could
>  transform to the needed Cartesian matrix?"
>
>   The input files I haved used in siesta calculation is attached as well
> as
> the
>output force-constant *.FC file.
>   We have used the executables fcbuild, siesta and vibrator included
> in
> the Util directory of
>   siesta-3.1 software package.First, we run
>$ fcbuild < graphene.fdf
>this builds a file FC.fdf which is a partial input for Siesta. Then, We
> run
>$ siesta < graphene-siesta.fdf
>Siesta provides the file graphene.FC with the force constant
> matrix.
>
>
>  Please help me.
>  best regards
>  bedamani
>

Re: [SIESTA-L] pdos of pz orbital

[apostnik]

>
> Hi dear siesta users,
>
> I plot pdos of pz orbital of carbon for benzene molecule which is
> functionalized with methyl group,
>
> in this curve I knew which energies is dominated by pz orbitals, but I
> want to know which energy is approximately related to pz orbital of
> methyle group. Does pdos help me?
>
> I really  need which energy appropriated to each atom in my typical
> system.
> Do I analyze incorrectly?

Dear Mehrzad,

plotting and looking at local pdos of atoms that interest you,
i.e. in the methyl group, may give some clue. However, it is likely
that they occupy certain interval of energies, overlapping with pdos
of other atoms.
In general, your question "which energy appropriated to each atom"
is not quite correctly posed. Each energy may correspond to a given
molecular orbital, most probably spread over several atoms.
The same atoms may enter differently into different MOs.
What you may wish to do is to visualize these orbitals, look at them
one by one and decide which one is mostly related to "pz of methyl group"
(plus probably something else), or whatever interests you.
For this, you'll calculate and plot wavefunctions. You have to decide
which functions (numbers from... to) are potentially interesting for you.
This idea you'll get from looking at the pdos, as I write at the beginning,
and comparing the energies with the list of eigenvalues.

Good luck

Andrei Postnikov

Re: [SIESTA-L] about O(N) method

[apostnik]

> Hi, siesta users:
> I noticed that one line in Siesta manual, saying:
> "The Ordern(N) subsystem is quite fragile and only
> works for systems with clearly separated occupied
> and empty states."
>
> Is there any obvious reason that it cannot handle
> metallic systems?
>
> Best regards,
> Bing
>

Sure, Bing, it is.
In metal, the states which are close to the Fermi energy may shift
back and forth between occupied and unoccupied domains. On each iteration,
the eigenvalues are ordered by their increased values, and the
AUFBAU PRINCIPLE applied -
the so many lowest states (counted over the Brillouin zone)
are declared OCCUPIED
(and contribute to the density in the next iteration),
and the rest thrown out.

The Order(N) approach (at list as it is implemented in Siesta, or more
general, as far as I remember) avoids diagonalization; instead it works
directly on orbitals or  density matrix elements, iterating them according
to a certain criterion. Hence the method does not know eigenvalues and
goes ahead for density, total energy, and other related properties.
(But does NOT provide the density of states, say).
Therefore the method cannot systematically order the eigenvalues
and find out where the occupied ones end. Instead, it needs to know this
in advance. In practice, in addition to the number of electrons
(which is known anyway), the chemical potential has to be provided
IN ADVANCE as a measure to identify the occupied states.
If you have a gap and safely put the chemical potential inside,
the calculation may go stable. If you have a metal, it ends up in a mess;
the iteration algorithms go astray.

Best regards

Andrei Postnikov

Re: [SIESTA-L] about O(N) method

[apostnik]

> Hi, siesta users:
> I noticed that one line in Siesta manual, saying:
> "The Ordern(N) subsystem is quite fragile and only
> works for systems with clearly separated occupied
> and empty states."
>
> Is there any obvious reason that it cannot handle
> metallic systems?
>
> Best regards,
> Bing
>

Sure, Bing, it is.
In metal, the states which are close to the Fermi energy may shift
back and forth between occupied and unoccupied domains. On each iteration,
the eigenvalues are ordered by their increased values, and the
AUFBAU PRINCIPLE applied -
the so many lowest states (counted over the Brillouin zone)
are declared OCCUPIED
(and contribute to the density in the next iteration),
and the rest thrown out.

The Order(N) approach (at list as it is implemented in Siesta, or more
general, as far as I remember) avoids diagonalization; instead it works
directly on orbitals or  density matrix elements, iterating them according
to a certain criterion. Hence the method does not know eigenvalues and
goes ahead for density, total energy, and other related properties.
(But does NOT provide the density of states, say).
Therefore the method cannot systematically order the eigenvalues
and find out where the occupied ones end. Instead, it needs to know this
in advance. In practice, in addition to the number of electrons
(which is known anyway), the chemical potential has to be provided
IN ADVANCE as a measure to identify the occupied states.
If you have a gap and safely put the chemical potential inside,
the calculation may go stable. If you have a metal, it ends up in a mess;
the iteration algorithms go astray.

Best regards

Andrei Postnikov

Re: [SIESTA-L] Phonon density of states

[apostnik]

Dear Bedamani

my simple code was written having in mind (big) supercells
where a fair approximation to the density of modes can be gained
from looking just at q=0 modes. No effort was included
to provide any integration over the Brillouin zone for the cases
q.ne.0. In order to prevent using it on such cases which would
lead to wrong results, the check you refer to was introduced.

I don't have a code which would exactly serve your purposes.
As a substitution, I have a script which expands the "partial"
FC file (i.e., that of a supercell with only "non-equivalent"
atoms shifted), into a full FC file (i.e., as if all atoms of supercell
have been shifted). In your case, it would mean giving a possibility
to calculate q=0 phonon in a 54-at. supercell without moving
each of 54 atoms, but only 2 inequivalent ones.

This script seems to work but was not yet thoroughly tested.

I don't think, however, that the density of modes of pure Si
is of enough special interest to be calculated this way.

Best regards

Andrei Postnikov



>  In order to calculate the frequency and phonon mode  here  i have used
> vibra
> package  which contains  two programs (fcbuild.f and vibra.f).
>
>
>
> I compiled this two  by using the command “make”.
>
>
>
> I used si54.fdf file as input (taken from /Util/vibra/examples  folder).
>
>
>
>-   in order to run fcbuild  i have used the command
>
> ./fcbuild < si54.fdf
>
>
>
>   this provide  a file FC.fdf which is a partial input for Siesta
>
>
>
>- After that I used input si54-siesta.fdf  (taken from
>/Util/vibra/examples  folder) and used the commamd:
>
>
>
> ./siesta < si54-siesta.fdf
>
>
>
> this provide the file si54.FC
>
>
>
>-  after that I  run vibrator  using the command
>
> ./vibrator < si54.fdf
>
>
>
>and get the
>
>-  output written in si54.bands which gives the  phonon band structure.
>- This also provide the si54.vectors file containing frequency and
> eigen
>vectors.
>
>
>
>
>
> Now in order to calculate the phonon density of state (ph dos)  i have
> used
> the code given by  - Dr. Andrei Postnikov – .  I have found it from
> internet
> where he has given it as a reply of a  mail of Sharat Chandra.
>  Here I have attach this in this folder  as ‘ph.f'.
>
> I compiled this by using the command 'gfortran ph.f -o phonon'
> To run it i have used the command
>   >./phonon
>
> then it ask me to Specify SystemLabel of .XV file, i have specified it (I
> just typed si54) .
> Then it ask me to Specify SystemLabel of .vectors file, i have specified
> it
> (I just typed si54).
>
> After specifying  the vectors file, it provides an error note. Which is
> given below
> 
>
>  '”Attention: q.ne.0 in .vectors file!
>
>   Error reading Eigenvector line, last iev=
> *
> Please help me to solve this problem.
>
> Thank you all
>

Re: [SIESTA-L] Wrong magnetic moement is coming for Mn-c6h6 chains

[apostnik]

> Thanks  Andre
>
> The structure is correct.

Good. I see that your benzene molecule is not much magnetic,
that is good as well.

> I have tried the same structure with Quantum
> Espresso and VASP. I am getting the correct magnetic moment of 1 uB/unit
> cell.

A priori I'd say the magnetic moment of 4.5 m_B at Mn atom
is not necessarily wrong, but it indicates a very small coupling
of Mn with the neighboring organics. If not physical,
this could be due to a too short basis –
which is minimal (SZ) in your case, anyway.
Take a DZP, as advised, with a not too large PAO energy shift.

As you have a different calculation at hand which you trust,
I'd in your place compare the details – where the bands / peaks
in the density of states are, what are their widths
(as there must be a dispersion along Z).
Do not concentrate just at the magnetic moment.

Best regards

Andrei Postnikov


> I am attachign the input and output files and the pseudopotentials.
> I have also tried the pseudopotentials given on siesta's website.
> But nothing changes.
>
> Please have a look on it.
>
> On Sun, Sep 2, 2012 at 12:33 AM,  wrote:
>
>> > Hi Siesta Users,
>> >
>> > I am trying to reproduce the results on Mn-C6H6 chains. But I am
>> > getting the magnetic moment 3 uB/unit cell rather than 1 uB/unit cell.
>> > I tried various pseudopotentails. But result is coming same.
>>
>> Well...
>>
>> > Can anybody help me.
>>
>> No,
>> unless you explain
>> - what you are doing, what you results (electronic structure etc.) are
>> and what do you expect instead (for me 3 m_B per Mn atom in an organic
>> a priori looks fine).
>> If a problem, it can be
>> - an error in structure data;
>> - an error in pseudoptential/basis (how about Mn3p? too short basis?);
>> - not an error at all, but different ways to define local magnetic
>> moment;
>> - an attempt to reproduce somebody's wrong results, or those obtained
>> with different XC potentiel, etc.
>>
>> Have a nice day
>>
>> Andrei Postnikov
>>
>>
>> > --
>> > Pankaj Kumar
>> > Ph. D. Scholar
>> > School of Basic Sciences
>> > Department of Physics
>> > IIT Mandi
>> > Mob. No. +91 9736694726
>> >
>>
>>
>
>
> --
> Pankaj Kumar
> Ph. D. Scholar
> School of Basic Sciences
> Department of Physics
> IIT Mandi
> Mob. No. +91 9736694726
>

Re: [SIESTA-L] Wrong magnetic moement is coming for Mn-c6h6 chains

[apostnik]

> Thanks  Andre
>
> The structure is correct.

Good. I see that your benzene molecule is not much magnetic,
that is good as well.

> I have tried the same structure with Quantum
> Espresso and VASP. I am getting the correct magnetic moment of 1 uB/unit
> cell.

A priori I'd say the magnetic moment of 4.5 m_B at Mn atom
is not necessarily wrong, but it indicates a very small coupling
of Mn with the neighboring organics. If not physical,
this could be due to a too short basis –
which is minimal (SZ) in your case, anyway.
Take a DZP, as advised, with a not too large PAO energy shift.

As you have a different calculation at hand which you trust,
I'd in your place compare the details – where the bands / peaks
in the density of states are, what are their widths
(as there must be a dispersion along Z).
Do not concentrate just at the magnetic moment.

Best regards

Andrei Postnikov


> I am attachign the input and output files and the pseudopotentials.
> I have also tried the pseudopotentials given on siesta's website.
> But nothing changes.
>
> Please have a look on it.
>
> On Sun, Sep 2, 2012 at 12:33 AM,  wrote:
>
>> > Hi Siesta Users,
>> >
>> > I am trying to reproduce the results on Mn-C6H6 chains. But I am
>> > getting the magnetic moment 3 uB/unit cell rather than 1 uB/unit cell.
>> > I tried various pseudopotentails. But result is coming same.
>>
>> Well...
>>
>> > Can anybody help me.
>>
>> No,
>> unless you explain
>> - what you are doing, what you results (electronic structure etc.) are
>> and what do you expect instead (for me 3 m_B per Mn atom in an organic
>> a priori looks fine).
>> If a problem, it can be
>> - an error in structure data;
>> - an error in pseudoptential/basis (how about Mn3p? too short basis?);
>> - not an error at all, but different ways to define local magnetic
>> moment;
>> - an attempt to reproduce somebody's wrong results, or those obtained
>> with different XC potentiel, etc.
>>
>> Have a nice day
>>
>> Andrei Postnikov
>>
>>
>> > --
>> > Pankaj Kumar
>> > Ph. D. Scholar
>> > School of Basic Sciences
>> > Department of Physics
>> > IIT Mandi
>> > Mob. No. +91 9736694726
>> >
>>
>>
>
>
> --
> Pankaj Kumar
> Ph. D. Scholar
> School of Basic Sciences
> Department of Physics
> IIT Mandi
> Mob. No. +91 9736694726
>

Re: [SIESTA-L] Wrong magnetic moement is coming for Mn-c6h6 chains

[apostnik]

> Hi Siesta Users,
>
> I am trying to reproduce the results on Mn-C6H6 chains. But I am
> getting the magnetic moment 3 uB/unit cell rather than 1 uB/unit cell.
> I tried various pseudopotentails. But result is coming same.

Well...

> Can anybody help me.

No,
unless you explain
- what you are doing, what you results (electronic structure etc.) are
and what do you expect instead (for me 3 m_B per Mn atom in an organic
a priori looks fine).
If a problem, it can be
- an error in structure data;
- an error in pseudoptential/basis (how about Mn3p? too short basis?);
- not an error at all, but different ways to define local magnetic moment;
- an attempt to reproduce somebody's wrong results, or those obtained
with different XC potentiel, etc.

Have a nice day

Andrei Postnikov


> --
> Pankaj Kumar
> Ph. D. Scholar
> School of Basic Sciences
> Department of Physics
> IIT Mandi
> Mob. No. +91 9736694726
>

Re: [SIESTA-L] Wrong magnetic moement is coming for Mn-c6h6 chains

[apostnik]

> Hi Siesta Users,
>
> I am trying to reproduce the results on Mn-C6H6 chains. But I am
> getting the magnetic moment 3 uB/unit cell rather than 1 uB/unit cell.
> I tried various pseudopotentails. But result is coming same.

Well...

> Can anybody help me.

No,
unless you explain
- what you are doing, what you results (electronic structure etc.) are
and what do you expect instead (for me 3 m_B per Mn atom in an organic
a priori looks fine).
If a problem, it can be
- an error in structure data;
- an error in pseudoptential/basis (how about Mn3p? too short basis?);
- not an error at all, but different ways to define local magnetic moment;
- an attempt to reproduce somebody's wrong results, or those obtained
with different XC potentiel, etc.

Have a nice day

Andrei Postnikov


> --
> Pankaj Kumar
> Ph. D. Scholar
> School of Basic Sciences
> Department of Physics
> IIT Mandi
> Mob. No. +91 9736694726
>

Re: [SIESTA-L] Wrong magnetic moement is coming for Mn-c6h6 chains

[apostnik]

> Hi Siesta Users,
>
> I am trying to reproduce the results on Mn-C6H6 chains. But I am
> getting the magnetic moment 3 uB/unit cell rather than 1 uB/unit cell.
> I tried various pseudopotentails. But result is coming same.
>
> Can anybody help me.

Sorry, I did not noticed that an input file
was attached to your mail.
An output would have been much more useful.
Still:
(in addition to my previous comments):
check that you initialize magnetism on Mn only
and NOT on C,H atoms (use DM.InitSpin block).
I guess you may get some magnetism where it does not belong.

Best regards

Andrei Postnikov


>
> --
> Pankaj Kumar
> Ph. D. Scholar
> School of Basic Sciences
> Department of Physics
> IIT Mandi
> Mob. No. +91 9736694726
>

Re: [SIESTA-L] Wrong magnetic moement is coming for Mn-c6h6 chains

[apostnik]

> Hi Siesta Users,
>
> I am trying to reproduce the results on Mn-C6H6 chains. But I am
> getting the magnetic moment 3 uB/unit cell rather than 1 uB/unit cell.
> I tried various pseudopotentails. But result is coming same.
>
> Can anybody help me.

Sorry, I did not noticed that an input file
was attached to your mail.
An output would have been much more useful.
Still:
(in addition to my previous comments):
check that you initialize magnetism on Mn only
and NOT on C,H atoms (use DM.InitSpin block).
I guess you may get some magnetism where it does not belong.

Best regards

Andrei Postnikov


>
> --
> Pankaj Kumar
> Ph. D. Scholar
> School of Basic Sciences
> Department of Physics
> IIT Mandi
> Mob. No. +91 9736694726
>

Re: [SIESTA-L] how to plot spin density

[apostnik]

Dear Yuly:

it is difficult to discuss without seeing your data.
Send me your .XSF file in case of doubt. But, in its essence
the matter is very simple. You have two blocks of data
(for spin-up and spin-down density), and you can explore them
with XCrySDen at your will. First, set the weight 1.0 for first block
and 0. for the second one. Is the result meaningful?
Then make it inversely, the weights 0. and 1.0. Is the spin-down
density meaningful? (It must be close to that for spin-up,
with some small deviations). If this is indeed the case
(none of the densities is zero everywhere, both are positive everywhere,
no crazy numbers), then the difference (weights 1.0 and -1.0)
MUST give you a reasonable spin density. You should keep in mind
that the magnitude of spin density is MUCH SMALLER than that of
charge density, so you should choose the isosurfaces correspondingly.
XCrySDen will show you the max. and min. numbers.
Sorry, I cannot much to that, without knowing what you see
and what you expect.

Best regards

Andrei Postnikov


> Dear Dr. Postnikov.
>
> Sorry I have made an mistake, What I mean that:
> When I choose subblock # 0 = 1.0 and subblock # 1 = -1.0, I got the result
> same as when I choose subblock # 0 = 1.0 and subblock # 1 = 1.0.
>
> The result is the plot of charge density like your picture in page 28 in
> your slide (visualization and post-processing tools for siesta).
>
> Previously, I hope that I will got the profil of spin density when I
> choose subblock # 0 = 1.0 and subblock # 1 = -1.0, as you said in the
> slide, page 26. "Choose 1.0 and 1.0 for plotting the charge density,1.0
> and −1.0
> for plotting the spin density".
>
> I am new in computational field. Please give me a clue.
>
>
> Thank you
>
> Best Regards
>
> Yuly Kusumawati
>
>
>
>
>>> Dear Siesta Users
>>>
>>> I have a problem when I have tried to plot spin density using XCrysden.
>>> I
>>> have followed the steps in the Dr. Postnikov slide's. And I have got
>>> the
>>> plot of charge density.
>>> When I tried to plot the spin density by change the multiply factor in
>>> subblock #1 = -0.1. What I have got is again the charge density, it is
>>> not
>>> the spin density that I want to plot.
>>
>> Dear Yuly:
>> I don't know what you expect by selecting
>>> subblock #1 = -0.1.
>> This gives you 10% of the spin-up density.
>> How about the spin-down one (subblock #2)?
>>
>> In order to get (spin-up) MINUS (spin-down),
>> take weights
>> +1.0 for subblock #1
>> AND
>> -1.0 for subblock #2
>> as, you say, was in the lecture.
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>>
>>> My questions are:
>>> 1. Is my step is right? Because I read on the Dr. Postnikov slide's in
>>> the
>>> page 12, that for plotting the spin density we have to choose 1.0 and
>>> -1.0
>>> in the subblock.
>>>
>>> 2. What is the right step for plotting spin density?
>>>
>>> Thank you for your answer.
>>>
>>>
>>> Best Regards
>>>
>>>
>>> Yuly Kusumawati
>>>
>>>
>>> --
>>>
>>> http://www.its.ac.id
>>>
>>>
>>
>>
>
>
>
>

Re: [SIESTA-L] how to plot spin density

[apostnik]

Dear Yuly:

it is difficult to discuss without seeing your data.
Send me your .XSF file in case of doubt. But, in its essence
the matter is very simple. You have two blocks of data
(for spin-up and spin-down density), and you can explore them
with XCrySDen at your will. First, set the weight 1.0 for first block
and 0. for the second one. Is the result meaningful?
Then make it inversely, the weights 0. and 1.0. Is the spin-down
density meaningful? (It must be close to that for spin-up,
with some small deviations). If this is indeed the case
(none of the densities is zero everywhere, both are positive everywhere,
no crazy numbers), then the difference (weights 1.0 and -1.0)
MUST give you a reasonable spin density. You should keep in mind
that the magnitude of spin density is MUCH SMALLER than that of
charge density, so you should choose the isosurfaces correspondingly.
XCrySDen will show you the max. and min. numbers.
Sorry, I cannot much to that, without knowing what you see
and what you expect.

Best regards

Andrei Postnikov


> Dear Dr. Postnikov.
>
> Sorry I have made an mistake, What I mean that:
> When I choose subblock # 0 = 1.0 and subblock # 1 = -1.0, I got the result
> same as when I choose subblock # 0 = 1.0 and subblock # 1 = 1.0.
>
> The result is the plot of charge density like your picture in page 28 in
> your slide (visualization and post-processing tools for siesta).
>
> Previously, I hope that I will got the profil of spin density when I
> choose subblock # 0 = 1.0 and subblock # 1 = -1.0, as you said in the
> slide, page 26. "Choose 1.0 and 1.0 for plotting the charge density,1.0
> and −1.0
> for plotting the spin density".
>
> I am new in computational field. Please give me a clue.
>
>
> Thank you
>
> Best Regards
>
> Yuly Kusumawati
>
>
>
>
>>> Dear Siesta Users
>>>
>>> I have a problem when I have tried to plot spin density using XCrysden.
>>> I
>>> have followed the steps in the Dr. Postnikov slide's. And I have got
>>> the
>>> plot of charge density.
>>> When I tried to plot the spin density by change the multiply factor in
>>> subblock #1 = -0.1. What I have got is again the charge density, it is
>>> not
>>> the spin density that I want to plot.
>>
>> Dear Yuly:
>> I don't know what you expect by selecting
>>> subblock #1 = -0.1.
>> This gives you 10% of the spin-up density.
>> How about the spin-down one (subblock #2)?
>>
>> In order to get (spin-up) MINUS (spin-down),
>> take weights
>> +1.0 for subblock #1
>> AND
>> -1.0 for subblock #2
>> as, you say, was in the lecture.
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>>
>>> My questions are:
>>> 1. Is my step is right? Because I read on the Dr. Postnikov slide's in
>>> the
>>> page 12, that for plotting the spin density we have to choose 1.0 and
>>> -1.0
>>> in the subblock.
>>>
>>> 2. What is the right step for plotting spin density?
>>>
>>> Thank you for your answer.
>>>
>>>
>>> Best Regards
>>>
>>>
>>> Yuly Kusumawati
>>>
>>>
>>> --
>>>
>>> http://www.its.ac.id
>>>
>>>
>>
>>
>
>
>
>

Re: [SIESTA-L] how to plot spin density

[apostnik]

> Dear Siesta Users
>
> I have a problem when I have tried to plot spin density using XCrysden. I
> have followed the steps in the Dr. Postnikov slide's. And I have got the
> plot of charge density.
> When I tried to plot the spin density by change the multiply factor in
> subblock #1 = -0.1. What I have got is again the charge density, it is not
> the spin density that I want to plot.

Dear Yuly:
I don't know what you expect by selecting
> subblock #1 = -0.1.
This gives you 10% of the spin-up density.
How about the spin-down one (subblock #2)?

In order to get (spin-up) MINUS (spin-down),
take weights
+1.0 for subblock #1
AND
-1.0 for subblock #2
as, you say, was in the lecture.

Best regards

Andrei Postnikov

>
> My questions are:
> 1. Is my step is right? Because I read on the Dr. Postnikov slide's in the
> page 12, that for plotting the spin density we have to choose 1.0 and -1.0
> in the subblock.
>
> 2. What is the right step for plotting spin density?
>
> Thank you for your answer.
>
>
> Best Regards
>
>
> Yuly Kusumawati
>
>
> --
>
> http://www.its.ac.id
>
>

Re: [SIESTA-L] how to plot spin density

[apostnik]

> Dear Siesta Users
>
> I have a problem when I have tried to plot spin density using XCrysden. I
> have followed the steps in the Dr. Postnikov slide's. And I have got the
> plot of charge density.
> When I tried to plot the spin density by change the multiply factor in
> subblock #1 = -0.1. What I have got is again the charge density, it is not
> the spin density that I want to plot.

Dear Yuly:
I don't know what you expect by selecting
> subblock #1 = -0.1.
This gives you 10% of the spin-up density.
How about the spin-down one (subblock #2)?

In order to get (spin-up) MINUS (spin-down),
take weights
+1.0 for subblock #1
AND
-1.0 for subblock #2
as, you say, was in the lecture.

Best regards

Andrei Postnikov

>
> My questions are:
> 1. Is my step is right? Because I read on the Dr. Postnikov slide's in the
> page 12, that for plotting the spin density we have to choose 1.0 and -1.0
> in the subblock.
>
> 2. What is the right step for plotting spin density?
>
> Thank you for your answer.
>
>
> Best Regards
>
>
> Yuly Kusumawati
>
>
> --
>
> http://www.its.ac.id
>
>

Re: [SIESTA-L] Phonon DOS

[apostnik]

Dear Mic,

without you having provided anything of your files,
I can only recommend you to look into your
.vectors file -
whether it is not empty, contains the data for 6 modes, etc.

Best regards

Andrei Postnikov


> Hello
> When I try Vibra(
> siesta-3.0-rc2/Util/Vibra/Examples/Si2
>  )to calculated the phDOS of Si using the phdos utility, I get the
> following message.
> How to fix this error.
>  
> Specify SystemLabel of .XV file: si2
>     1 different types found in the XV file:
>     2 atoms of type 1
>  Specify SystemLabel of .vectors file: si2
>   Error reading Eigenvector line, last iev=  
>  
> Mic
>
>
> --- On Wed, 5/31/06, Andrei Postnikov  wrote:
>
>
> From: Andrei Postnikov 
> Subject: Re: [SIESTA-L] Phonon DOS
> To: siest...@listserv.uam.es
> Date: Wednesday, May 31, 2006, 2:46 AM
>
>
> On Wed, 31 May 2006, Sharat Chandra wrote:
>
> | Hi
> | Can any one please tell me how to calculate the phonon DOS using Vibra
> | module? Has any one developed a code for carrying out the calculation?
>
> Hallo Sharat:
> yes I have a tool to calculate (local) phonon DOS,
> however it is of limited usefulness:
> it assumes a single q-point in a presumable large supercell,
> and smears out the peaks at frequencies (and with weights)
> as provided by the "vibrator".
> Please find enclosed the source code and README.
> Best regards,
>
> Andrei
>
> +-- Dr. Andrei Postnikov  Tel. +33-387315873 - mobile
> +33-666784053 ---+
> | Paul Verlaine University - Institute de Physique Electronique et
> Chimie,     |
> | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz,
> France |
> +-- apost...@uos.de 
> http://www.home.uni-osnabrueck.de/apostnik/ --+

Re: [SIESTA-L] Phonons negative frequencies

[apostnik]

> Andrei,
>
> thank you very much for your answer. I had already found in the ML
> database
> the nice phonon reference you indicated to me.
> By negative frequencies, I mean:
>
>   eigenvalue #   1  omega=  -40.738952961112830
>   eigenvalue #   2  omega= -2.05736482791876937E-002
>   eigenvalue #   3  omega= -1.32180970904245670E-002

Dear Carlo,
Your frequencies 2 and 3 are zero, for all practical use.
No special importance therefore that they are negative.
If you have yet the third frequency (positive, in your case?)
which is as close to zero as those two (your Nr4?)
then likely you don't have a problem with MeshCutoff.
However, your frequency Nr1 is an indication
that your relaxation is not good enough.
The negative frequency (if a stable one and not due to numerical noise)
means: as you displace the atoms in some combination, the energy
goes down. It shouldn't be like this at the equilibrium.
Have a look into p.18 of
http://www.home.uni-osnabrueck.de/apostnik/Lectures/APostnikov-Phonons.pdf

> My system is a molecular junction. Just to explain you, at the beginning I
> calculated the electrodes bulk lattice parameter, and then I built the
> junction geometry. If I then relax by setting MD.variableCell=T the cell
> dimensions increase by 1.5%... Do you think that I should necessarily use
> a variable cell in the relaxation? (I am suspicious about it)

Not necessarily at all; you can calculate phonons
in a constrained system, under pressure etc. However,
atom positions have to be optimized subject to this constraint
(to be sure that trial atom displacements are around the equilibrium).

> Finally, could you please give me some reference about the eggbox effect
> test?

It has been discussed many many times in the main list and in tutorials.
E.g., you may have a look in
http://www.home.uni-osnabrueck.de/apostnik/Lectures/SIESTA-tuto.pdf
(pp.10-12)

Good luck

Andrei Postnikov

>
> Thank you very much!
>
> Carlo
>
>
> 2012/4/7 
>
>> Dear Carlo,
>> as you don't specify how negative your negative frequencies are,
>> it is not so easy to judge... Your sum of forces seems a bit too high...
>> Even if not dramatically... And some individual forces as well -
>> e.g., Fx on atom 2. So generally I'd recommend even higher MeshCutoff
>> and further relaxation. A priori 350 Ry is not safely high,
>> sometimes one needs to go to much higher.
>>
>> Before doing the phonons, I'd recommend the usual "eggbox test"
>> to be sure that fluctuations of summary force as function of
>> uniform displacement of all atoms remain within |0.1|.
>>
>> Concerning specifically the phonons, you may wish to check
>> http://www.home.uni-osnabrueck.de/apostnik/Lectures/APostnikov-Phonons.pdf
>> - especially p.18 of it.
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>> > Dear Siesta users,
>> >
>> > I am trying to calculate phonon frequencies at Gamma for a molecular
>> > junction.
>> > I relaxed the system by keeping fixed the Lattice Vectors (not a
>> variable
>> > cell). When I compute the frequencies (with MeshCutoff = 200 Ry), I
>> obtain
>> > that the first three frequencies are negative, which indicates that
>> the
>> > results are not correct.
>> > Then, if I increase the MeshCutoff to 350 Ry (which seems a high value
>> to
>> > me) I obtain again negative frequencies.
>> > Could you suggest me how to fix the problem?
>> > The following are the forces of the system...do you think that they
>> are
>> > too
>> > high? Which is a good criteria to say that the forces are low enough?
>> >
>> > Thank you very much,
>> >
>> > Carlo
>> >
>> > siesta:  10.106233   -0.0961800.047641
>> > siesta:  2   -0.6558380.034936   -0.004086
>> > siesta:  30.364736   -0.1003680.166114
>> > siesta:  40.1306900.010500   -0.101584
>> > siesta:  5   -0.1224100.420815   -0.012044
>> > siesta:  60.120625   -0.0312920.152270
>> > siesta:  70.0033190.004051   -0.012274
>> > siesta:  80.004951   -0.0157190.001451
>> > siesta:  9   -0.0095230.0046100.013162
>> > siesta: 100.019411   -0.008606   -0.009840
>> > siesta: 110.002485   -0.0101310.005577
>> > siesta: 120.001815   -0.012882   -0.010601
>> > siesta: 130.002293   -0.0278280.004317
>> > siesta: 140.003541   -0.004099   -0.003595
>> > siesta: 15   -0.004545   -0.005597

Re: [SIESTA-L] siesta spin polarization problem

[apostnik]

>  Hi,all
> I use the siesta to calculate the spin polarizaiton of a Ti doped
> system without magnetic materials. When no spin-polarized, the total
> energy is  -9276.290371 eV; while adding the spin polarizaiton, the
> total energy is -9279.474363 eV, and it shows "Total spin polarization
> (Qup-Qdown) =4.00" in the output file. My question are

> :1. Do
> the calcaluted total energy trend is right, adding the spin
> polarization causing the absolute value bigger?

Dear Shi:
The total energy goes down with magnetism, this seems reasonable
(if the magnetic solution, indeed, is the stable one)

>2. Is it possible that
> Qup-Qdown value is "4"?
Everything is possible. It is not quite clear what you are calculating,
"a Ti doped system without magnetic materials",
but you should check whether the magnetism comes from where it is
expected to: from 2 d-electrons on each Ti atom (do you have two of them
in your system?) and not from something else strangely polarized.

>3. Which unit is the "Total spin polarization"?
Bohr magnetons, i.e., difference in (spin-up)-(spin-down)
electron numbers.

Best regards

Andrei Postnikov

>   
>  Shi
>   
>  2012.4.6
>

Re: [SIESTA-L] Phonons negative frequencies

[apostnik]

Dear Carlo,
as you don't specify how negative your negative frequencies are,
it is not so easy to judge... Your sum of forces seems a bit too high...
Even if not dramatically... And some individual forces as well -
e.g., Fx on atom 2. So generally I'd recommend even higher MeshCutoff
and further relaxation. A priori 350 Ry is not safely high,
sometimes one needs to go to much higher.

Before doing the phonons, I'd recommend the usual "eggbox test"
to be sure that fluctuations of summary force as function of
uniform displacement of all atoms remain within |0.1|.

Concerning specifically the phonons, you may wish to check
http://www.home.uni-osnabrueck.de/apostnik/Lectures/APostnikov-Phonons.pdf
- especially p.18 of it.

Best regards

Andrei Postnikov


> Dear Siesta users,
>
> I am trying to calculate phonon frequencies at Gamma for a molecular
> junction.
> I relaxed the system by keeping fixed the Lattice Vectors (not a variable
> cell). When I compute the frequencies (with MeshCutoff = 200 Ry), I obtain
> that the first three frequencies are negative, which indicates that the
> results are not correct.
> Then, if I increase the MeshCutoff to 350 Ry (which seems a high value to
> me) I obtain again negative frequencies.
> Could you suggest me how to fix the problem?
> The following are the forces of the system...do you think that they are
> too
> high? Which is a good criteria to say that the forces are low enough?
>
> Thank you very much,
>
> Carlo
>
> siesta:  10.106233   -0.0961800.047641
> siesta:  2   -0.6558380.034936   -0.004086
> siesta:  30.364736   -0.1003680.166114
> siesta:  40.1306900.010500   -0.101584
> siesta:  5   -0.1224100.420815   -0.012044
> siesta:  60.120625   -0.0312920.152270
> siesta:  70.0033190.004051   -0.012274
> siesta:  80.004951   -0.0157190.001451
> siesta:  9   -0.0095230.0046100.013162
> siesta: 100.019411   -0.008606   -0.009840
> siesta: 110.002485   -0.0101310.005577
> siesta: 120.001815   -0.012882   -0.010601
> siesta: 130.002293   -0.0278280.004317
> siesta: 140.003541   -0.004099   -0.003595
> siesta: 15   -0.004545   -0.0055970.018501
> siesta: 16   -0.0023750.005693   -0.016887
> siesta: 170.007410   -0.017539   -0.001118
> siesta: 180.015827   -0.012960   -0.005705
> siesta: 19   -0.010629   -0.007460   -0.001509
> siesta: 200.0087480.000683   -0.004698
> siesta: 21   -0.0049430.029703   -0.008064
> siesta: 22   -0.0032020.0045000.030359
> siesta: 23   -0.009982   -0.003884   -0.044988
> siesta: 240.012110   -0.008463   -0.054600
> siesta: 250.0148910.0001140.030375
> siesta: 260.0117140.007854   -0.002981
> siesta: 27   -0.006905   -0.0246510.005841
> siesta: 280.002399   -0.0193180.013282
> siesta: 29   -0.0198520.006114   -0.030653
> siesta: 30   -0.0236870.004007   -0.028044
> siesta: 310.023195   -0.022556   -0.014991
> siesta: 320.024333   -0.033430   -0.028258
> siesta: 33   -0.0149370.0063780.009843
> siesta: 340.004912   -0.019004   -0.008074
> siesta: 350.006259   -0.0123230.009612
> siesta: 36   -0.0104530.0172220.013853
> siesta: 370.0012810.001284   -0.002971
> siesta: 380.0095640.017869   -0.013170
> siesta: 390.004626   -0.016715   -0.008634
> siesta: 400.008447   -0.028871   -0.007756
> siesta: 41   -0.020042   -0.027721   -0.021694
> siesta: 420.020027   -0.019833   -0.095505
> siesta: 
> siesta:Tot0.016519   -0.011098   -0.032128
>
>
> --
>
> Carlo Motta
> * Ph.D. Student in Materials Science*
> University Office Phone : +39 02 6448-5183
> Department of Materials Science
> University of Milano-Bicocca
> Via R. Cozzi 53, 20125  Milano,  Italy
>

Re: [SIESTA-L] BZ sampling with user generated list of kpoints

[apostnik]

> Hello,
>
> Is there any option is siesta/transiesta to sample the BZ only at
> specific, user-specified kpoints? The manual gives only two options:
> Monkhorst Pack grid and kgrid cutoff.

Dear Diana -
these are the options provided for full BZ integration,
so it's logical that they are taken somehow regular, no?
Otherwise - if you do not make iterations
but want just sample data over selected k-points,
you can proceed exactly the same way as for band plotting:
you define the list of k-points at you wish
("band lines" scanning the vicinity of your favourite k-point)
and do with the results what you want...

Best results

Andrei Postnikov

> For transport calculations I am
> interested in sampling only a small area around a high symmetry point,
> where I expect most of the transport to take place.
>
> Thank you for your help,
>
> Diana
>
>
> Diana Otálvaro
> d.otalv...@utwente.nl
> Computational Material Science
> MESA+ Institute for Nanotechnology
> University of Twente
> Carre 4049
> Postbus 217
> NL-7500 AE Enschede
> tel: +31-53-489-2986

Re: [SIESTA-L] Band Structure Results

[apostnik]

>
>> And another things that I want to ask is, how to label the x-axis so
>> that
>> we get the symmetry label on that (ex. Gamma, M, K etc.)
>
> I think they are not automatically produced by the Gnubands
> script, but knowing where they are on the abscissa axis
> (that's you who defined the number of k-points between
> the symmetric corner points, right?),
> you can hopefully add them by hand in your favourite plotting routine.

Dear Yuly,
a correction to this:
I think, the last lines in the .bands file
give the direct placement of labels
along the abscissa axis of the plot.

So in the worst case you just add them by hand
at the specified points.

Best regards

Andrei Postnikov

Re: [SIESTA-L] Finding Magnetic Moments Using SIESTA

[apostnik]

> Hello,
> I have recently begun working with SIESTA to study the Iron Pncitide
> LaOFeAs.  I have successfully obtained a few converged SCF runs and have
> used DenChar to plot magnetisation (the difference between up and down
> spins in LSDA) through various planes in the unit cell.  However, rather
> than the values on these planes I am interested in the integrated moment
> around certain atoms within my unit cell so that I may compare my
> results with experimental values (c.f.
> http://arxiv.org/pdf/1001.0536.pdf).  Does anybody know of a way of
> obtaining such integrated values? Thank You
> Jack
>
> P.s. If anyone knows what dDmax means in the SCF iterations I have been
> searching but still cant work it out!
>


Dear Jack,
I hope my posting from exactly today 5 years back
http://www.mail-archive.com/siesta-l@listserv.uam.es/msg01232.html
may be helpful;
the tool in question is still there:
http://www.home.uni-osnabrueck.de/apostnik/Software/grdint.f
Please note that it makes the most primitive summation
over the grid points within a given sphere, and is more slow
than ought to be. An example of its use (comparison to VASP)
is in
DOI:10.1103/PhysRevB.82.054418
Enjoy / fee free to modify

Best regards

Andrei Postnikov

Re: [SIESTA-L] Band Structure Results

[apostnik]

> Hallo SIESTA users!
>
> I have run my fdf to get the band structure for graphene and TIO2 Rutile.
> I have made the visualization using gnubands and gnuplot. But then, I have
> very much lines. Can anyone informs me why it can be happened? is there
> any mistake?

Dear Yuly,
it is difficult to tell without knowing details about
your input/output, by looking just as your figures.
The first impression is, you have many bands
due to using of an enlarged supercell.
A priori, you can easily estimate
the expected number of occupied bands, simply counting
1*s, 3*p etc. for occupied states
which enter as valence (in your pseudopotential),
summing this up over all atoms of the unit cell.

> And another things that I want to ask is, how to label the x-axis so that
> we get the symmetry label on that (ex. Gamma, M, K etc.)

I think they are not automatically produced by the Gnubands
script, but knowing where they are on the abscissa axis
(that's you who defined the number of k-points between
the symmetric corner points, right?),
you can hopefully add them by hand in your favourite plotting routine.

Otherwise you can be inspired by more sophisticated scripts
from other band structure codes (WIEN2k, TB-LMTO, ...);
this part of work has nothing to do with the specifics of SIESTA.

Best regards

Andrei Postnikov


>
> I have attached the file of bands structure visualization
> and here the K-point for graphene
> %block Bandlines
> 1  -0.500   0.000  0.000 M
> 25  0.000   0.000  0.000 \Gamma
> 25  0.3329 -0.6658 0.000 K
> %endblock BandLines
>
> while for tio2 (rutile, space group 136)
> %block Bandlines
> 1  0.000 0.000 0.000 \Gamma
> 20 0.000 0.500 0.000 X
> 25 0.000 0.500 0.500 R
> 20 0.000 0.000 0.500 Z
> 25 0.000 0.000 0.000 \Gamma
> 20 0.500 0.500 0.000 M
> 25 0.500 0.500 0.500 A
> 20 0.000 0.000 0.500 Z
> %endblock BandLines
>
> Thank you,
>
> Best Regards
>
> Yuly Kusumawati
>
> Laboratory of Computational Chemistry
> Institut Teknologi Bandung (ITB)
>
>
> --
>
> http://www.its.ac.id
>

Re: [SIESTA-L] xv2xsf conversion

[apostnik]

> Dear Siesta Users...
>
> I want to convert the *.XV file to *.XSF file which is readable
> byXCrySDen.
> I manage to compile the codes in siesta-3.1/Util/Contrib/Apostnikov. when
> i
> convert the file using xv2xsf, i get the file in *.XSF format but i get
> the
> following error:
>
> Specify  SystemLabel (or 'siesta' if none):   Error opening file
> 0.0.XV as old formatted
>
> same i get in *.XSF file..
>
> is there anybody who is able to help me?

Yes, dear Pradeep.
you start the xv2xsf script without any parameters;
when prompted
"Specify  SystemLabel (or 'siesta' if none):"
you type your SystemLabel (the first part of the name
of your .XV file, that comes in from of ".XV")
and press "Enter"

Best regards

Andrei Postnikov



> thanks in advance
> --
> Pradeep Kumar
> India
>

Re: [SIESTA-L] phonon calculation error

[apostnik]

Dear Fenhong,
you are right,
"Siesta can only calculate frequency of whole system",
and that's for a good reason: vibration in a coupled system
is a property of the whole system; remember an exercise
on coupled oscillators in the first year in physics?
This does not mean that you cannot ask a question,
"at which frequency this particular bond vibrates"? -
It may vibrate at many frequencies, involving its other neighbors
to bigger or smaller extent. You go through analysis of
collective vibration modes, and you make your conclusions.
I think, you want to can ask a question differently:
"How this particular part of system will vibrate if decoupled
from the rest of molecule"? - This is something that is not
possible to probe in experiment, but in computer experiment
it is legitimate. Only that you have to clearly define
your physical model. I can imagine 3 possibilities:
1) The other atoms do not exist and react in no way
on the fragment of my interest. Then you throw these atoms away
and calculate just your small fragment (if it still makes sense...)
by Siesta and then Vibra. Most probably, this won't have much sense
because of wrong chemistry.
2) The other atoms exist in what regards the electronic structure,
but you want to neglect the forces which exist between them
and your selected fragment, in what regards the calculation
of vibrations. Then, you need to calculate FC only
displacing the atoms within your selected fragment, but then
you'd have to tweak the .FC file by hand, removing from it
the interactions with all external atoms. This is very easy
once you know how the .FC file is organized. Then you'd need
to tweak the Vibra input file, so that it would now relate
only to your reduced fragment, throwing out the rest.
(If you calculate a molecule and hence no dispersion is involved,
you can just manipulate the .XV file and use my simple "emulator"
of Vibra, but that's a different story).
Keep in mind that this is a very special and very unphysical
assumption, to switch these interactions away... But certainly
you good right to do, as long as you know what you are doing.
3) You take into account all neighboring atoms and interactions
with them, but you don't want these atoms to move. They'll remain
immobile and exercise forces on your selected atoms as they should,
from Siesta calculation. For this model, you calculate FC over
all atoms, but in the Vibra calculations, you assign infinite
(in practical terms, very big) masses to all atoms but your
selected ones. Then your selected fragment will vibrate like
in a cage.

Beyond these three possibilities, I don't see which physical model
you might have in mind.

My excuses for a long message.

Best regards

Andrei Postnikov



> Dear Siesta developers and users,
> Is it possible to only calculate part of system?
> I am going to calculate phonon frequency of   part of the molecule.
> There are 13 atoms in molecule. so I fix No 11 and 12 by
> "MD.FCfirst  11
> MD.FClast  12
> "
> When I use "Vibra" to calculate the vibration frequency between atom 11
> and
> 12,  It seems there is something wrong with the *.FC file
> here is the msg
> "lxmax= 0
> lymax= 0
> lzmax= 0
> Number of unit cells in Supercell  = 1
> Eigenvectors =   True
> Computing Eigenvalues and Eigenvectors
> forrtl: severe (24): end-of-file during read, unit 11, file H.FC"
> vibrator   0049780A  Unknown   Unknown
> Unknown
> vibrator   00496A0A  Unknown   Unknown
> Unknown
> vibrator   00457F52  Unknown   Unknown
> Unknown
> vibrator   0041F531  Unknown   Unknown
> Unknown
> vibrator   0041EE1F  Unknown   Unknown
> Unknown
> vibrator   0043E223  Unknown   Unknown
> Unknown
> vibrator   004050D2  Unknown   Unknown
> Unknown
> vibrator   00402F82  Unknown   Unknown
> Unknown
> libc.so.6  00396D21D8B4  Unknown   Unknown
> Unknown
> vibrator   00402EA9  Unknown   Unknown
> Unknown
>
>
> Does anyone know how to solve it or Siesta can only calculate frequency of
> whole system  ?
> Thanks in advance.
>
> Regard,
> Fenhong
>

Re: [SIESTA-L] Compress system using SIESTA??

[apostnik]

In principle, yes... But take care:
if you work over a large interval of pressures
with a fixed-basis code like Siesta, the quality of basis
is not equally good everywhere,
so that you might be pretty off either on the low-pressure end
or on the high-pressure end, depending on how/where you basis
has been optimized.
A test on a smaller representative system
against a planewave code might be useful,
to estimate the error...

Best regards

Andrei Postnikov


> Of course. Assuming your start denisity is not to faw away from a real
> system, you can use variable cell optimization, and increase the target
> pressure. If you have 2 values for the pressure, you could could estimate
> by linear extrapolation the final pressure which corresponds to the target
> density.
>
> Best regards
> Michael
>> hello!!
>>
>> i have a long polymer chain of about 420 atoms ..
>> can i compress it to target density of 1.35gm/cm3  using SIESTA??
>> As i have seen in a paper in which people take a long chain of polymer
>> and
>> wrap around cluster and then compress the system to target(bulk density)
>> using  GROMACS.
>>
>> Can it be possible with SIESTA? if yes then please tell me how?
>>
>> thanks ..
>>
>
>
> --
> Michael Schuch
> FB. Physik / AG. Statistische Physik
> Barbarastrasse 7
> 49069 Osnabrück
>
>

Re: [SIESTA-L] Platina is non-magnetic or what it the ATM_POP?

[apostnik]

Thanks Marcos,
you are right.
Then, it must be the dimer issue which kills the magnetism?
It would make sense to stop after some few iterations
and have a look at the magnetic moments
(and at the structure, e.g. the XV file, as well)

Best regards

Andrei


> Hi Andrei,
>
> Just one comment. If I'm not mistaken, the default spin polarization for
> Siesta, when you don't specify anything, is that one atom is initialized
> with spin up and the other down, in this case - an antiferro
> configuration.
> The final spin of the system could be non zero in this case, I suppose.
>
> Furthermore, I was checking something on Pt and Pd some days ago on the
> siesta mailing list, to see what people had done, and it seems that Pt is
> a
> bordeline case for magnetism, when it comes to wires. Results would differ
> with the use of GGA and LDA - check work by Anna Delin and Erio Tosatti on
> PRL, as well as the associated comment by Simone Alexandre and Jose Soler,
> and the reply to the comment by the authors of the paper.
>
> I just gave it a bird's eye look on the subject, and I recommend a more
> careful look at the siesta mailing list and the corresponding papers.
>
> Best regards,
>
> Marcos
>
> On Fri, Dec 2, 2011 at 8:44 AM,  wrote:
>
>> Dear isivkov,
>>
>> You use
>> SpinPolarized T#default
>> but I don't see if you ever made spin up different
>> from spin down (by InitSpin or otherwise).
>> By default, they start equal and remain equal.
>>
>> Moreover:
>> I don't see that you defined
>> AtomicCoordinatesFormat
>> so it must be Bohr by default? -
>> then you probably have Pt2 dimers at 1 Bohr distance,
>> separated by 20 Ang.
>>
>> In addition, i seems weird
>> to construct an empty box of 120 Ang size
>> (in X and Y dimensions), but this has nothing to do
>> with the magnetism issue.
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>> > Hello all.
>> >
>> > I decided to calculate a single atom of platinum, Not exactly single,
>> but
>> > chain of platinum with distance between atoms 10 A. I think it is like
>> a
>> > single atom. My .fdf file is here below. It is very strange, that
>> platinum
>> > is non-magnetic (as you can see from output file below), wile this is
>> the
>> > single atom. Atom of platinum must have magnetic moment.
>> >
>> > My input file
>> > ==
>> > #
>> >
>> -
>> > # FDF for Pt bulk
>> > #
>> > # LDA
>> > # Scalar-relativistic pseudopotential with non-linear partial-core
>> > correction
>> > #
>> > #
>> >
>> -
>> >
>> >  Cu bulk ###
>> >
>> > SystemName  PtLinChain
>> > SystemLabel PtLinChain   # Short name for naming files
>> >
>> > # Output options
>> >
>> > WriteCoorStep  true
>> >
>> > WriteMullikenPop   1
>> >
>> > # Species and atoms
>> >
>> > NumberOfSpecies1
>> > NumberOfAtoms  2
>> >
>> > %block ChemicalSpeciesLabel
>> >   1  78  Pt.LDA
>> > %endblock ChemicalSpeciesLabel
>> >
>> >
>> >
>> > %block PAO.Basis
>> >Pt.LDA 2 split 0.00  # Species label, number of l-shells
>> >n=6 0 2 P 1  # n, l, Nzeta, Polarization, NzetaPol
>> >0.00 0.00# 0.0 => default [6.982  5.935 \n 1.000  1.000]
>> >n=5 2 2  # n, l, zeta
>> >0.00 0.00
>> > %endblock PAO.Basis
>> >
>> >
>> >
>> > LatticeConstant  10  Ang
>> >
>> > %block LatticeVectors
>> >   12.000  0.000  0.000
>> >   0.000  12.000  0.000
>> >   0.000  0.000   2.000
>> > %endblock LatticeVectors
>> >
>> > %block kgrid_Monkhorst_Pack
>> >   1  0  0  0
>> >   0  1  0  0
>> >   0  0  8  0
>> > %endblock kgrid_Monkhorst_Pack
>> >
>> >
>> > XC.functional LDA   # Exchange-correlation functional
>> > XC.authorsCA# Exchange-correlation version
>> >
>> > MeshCutoff   150 Ry# Mesh cutoff. real space mesh
>> >
>> > # SCF options
>> > MaxSCFIterations  200   # Maximum number of SCF iter
>> > DM.MixingWeight   0.02   # New DM amount for next SCF
>> cycle
>> > DM.Tolerance  1.d-4 # Tolerance in maximum difference
>> > # between input and output DM
>> > DM.UseSaveDM  true  # to use continuation files
>> > DM.NumberPulay 5
>> >
>> > Diag.DivideAndConquer  .false.
>> > SolutionMethod diagon   # OrderN or Diagon
>> > ElectronicTemperature  25 meV   # Temp. for Fermi smearing
>> >
>> > SpinPolarized T#default
>> >
>> > # MD options
>> > #MD.TypeOfRun   cg   # Type of dynamics:
>> >
>> > #MD.UseSaveCG .true.  # for restarting
>> > #MD.UseSaveXV F   # atomic coords
>> >
>> > #MD.NumCGsteps  0# Number of CG steps for
>> > #   coordinate optimization
>> >

Re: [SIESTA-L] Platina is non-magnetic or what it the ATM_POP?

[apostnik]

Dear isivkov,

You use
SpinPolarized T#default
but I don't see if you ever made spin up different
from spin down (by InitSpin or otherwise).
By default, they start equal and remain equal.

Moreover:
I don't see that you defined
AtomicCoordinatesFormat
so it must be Bohr by default? -
then you probably have Pt2 dimers at 1 Bohr distance,
separated by 20 Ang.

In addition, i seems weird
to construct an empty box of 120 Ang size
(in X and Y dimensions), but this has nothing to do
with the magnetism issue.

Best regards

Andrei Postnikov


> Hello all.
>
> I decided to calculate a single atom of platinum, Not exactly single, but
> chain of platinum with distance between atoms 10 A. I think it is like a
> single atom. My .fdf file is here below. It is very strange, that platinum
> is non-magnetic (as you can see from output file below), wile this is the
> single atom. Atom of platinum must have magnetic moment.
>
> My input file
> ==
> #
> -
> # FDF for Pt bulk
> #
> # LDA
> # Scalar-relativistic pseudopotential with non-linear partial-core
> correction
> #
> #
> -
>
>  Cu bulk ###
>
> SystemName  PtLinChain
> SystemLabel PtLinChain   # Short name for naming files
>
> # Output options
>
> WriteCoorStep  true
>
> WriteMullikenPop   1
>
> # Species and atoms
>
> NumberOfSpecies1
> NumberOfAtoms  2
>
> %block ChemicalSpeciesLabel
>   1  78  Pt.LDA
> %endblock ChemicalSpeciesLabel
>
>
>
> %block PAO.Basis
>Pt.LDA 2 split 0.00  # Species label, number of l-shells
>n=6 0 2 P 1  # n, l, Nzeta, Polarization, NzetaPol
>0.00 0.00# 0.0 => default [6.982  5.935 \n 1.000  1.000]
>n=5 2 2  # n, l, zeta
>0.00 0.00
> %endblock PAO.Basis
>
>
>
> LatticeConstant  10  Ang
>
> %block LatticeVectors
>   12.000  0.000  0.000
>   0.000  12.000  0.000
>   0.000  0.000   2.000
> %endblock LatticeVectors
>
> %block kgrid_Monkhorst_Pack
>   1  0  0  0
>   0  1  0  0
>   0  0  8  0
> %endblock kgrid_Monkhorst_Pack
>
>
> XC.functional LDA   # Exchange-correlation functional
> XC.authorsCA# Exchange-correlation version
>
> MeshCutoff   150 Ry# Mesh cutoff. real space mesh
>
> # SCF options
> MaxSCFIterations  200   # Maximum number of SCF iter
> DM.MixingWeight   0.02   # New DM amount for next SCF cycle
> DM.Tolerance  1.d-4 # Tolerance in maximum difference
> # between input and output DM
> DM.UseSaveDM  true  # to use continuation files
> DM.NumberPulay 5
>
> Diag.DivideAndConquer  .false.
> SolutionMethod diagon   # OrderN or Diagon
> ElectronicTemperature  25 meV   # Temp. for Fermi smearing
>
> SpinPolarized T#default
>
> # MD options
> #MD.TypeOfRun   cg   # Type of dynamics:
>
> #MD.UseSaveCG .true.  # for restarting
> #MD.UseSaveXV F   # atomic coords
>
> #MD.NumCGsteps  0# Number of CG steps for
> #   coordinate optimization
> #MD.MaxCGDispl  0.05 Ang  # Maximum atomic displacement
> #   in one CG step (Bohr)
> #MD.MaxForceTol 0.005 eV/Ang  # Tolerance in the maximum
> #   atomic force (Ry/Bohr)
>
> # Atomic coordinates
> AtomicCoordinatesFormat ScaledCartesian
>
> # %block Zmatrix
> #cartesian
> #1 0.  0.  0.  1 1 1
> #1 0.3535  0.3535  0.5000  1 1 1
> #1 0.  0.  1.  1 1 1
> #1 0.3535  0.3535  1.5000  1 1 1
> #1 0.  0.  2.  1 1 1
> #1 0.3535  0.3535  2.5000  1 1 1
> #1 0.  0.  3.  0 0 0
> #1 0.3535  0.3535  3.5000  0 0 0
> #1 0.  0.  4.  0 0 0
> #1 0.3535  0.3535  4.5000  0 0 0
> # %endblock Zmatrix
>
> %block AtomicCoordinatesAndAtomicSpecies
> 0. 0. 0.000 1
> 0. 0. 1.000 1
> %endblock AtomicCoordinatesAndAtomicSpecies
>
> #%block GeometryConstraints
> #position from 1 to 4
> #%endblock GeometryConstraints
> =
>
>
> My output file(parts):
> At first SIESTA takes valence configuration from .psf file,
> as was in .inp file.
> 
> Reading pseudopotential information in formatted form from Pt.LDA.psf
>
> Pseudopotential generated from an atomic spin-polarized calculation
>
> Valence configuration for pseudopotential generation:
> 6s(1.00,0.00) rc: 2.32
> 6p(0.00,0.00) rc: 2.47
> 5d(5.00,4.00) rc: 1.23
> 5f(0.00,0.00) rc: 2.32
> For Pt.LDA, standard SIESTA heuristics set lmxkb to 3
>  (one more than the basis l, including polarization orbit

Re: [SIESTA-L] about the band structure calculation

[apostnik]

>
> Dear everyone,
> I try to calculate the band structure of ta2o5, but it is too
> time-consumption, so i doubt  that the input file is something wrong,
> please see in the attachment and help me check it.Thanks very much!
> Bo  Xiao

Dear Bo Xiao:
your sustem is big and periodic (lot of overlaps);
nobody knows how big your computer is.
Are you happy with your "regular" electronic structure calculation
(on a k-mesh of your ~130 points)?
If yes, then the band structure is nothing special,
it's just a calculation with different k-points (~300 in your case) ...
If it takes too long you can cut in in pieces;
you don't have to run it in a single run.

Beside the point: are you sure SZP is good for your purposes?

Best regards

Andrei Postnikov

Re:Re: [SIESTA-L] Negative value in phonon spectrum

[apostnik]

> Dear Postnikov,
>Thanks for your help ! As your advice, I changed the supercell (3 3
> 0) into (4 4 0), however, new problem occured. I tried many
> parameters, but there was always the same error:
>   "+  Segmentation fault "
> I don't konw what this error means and how to eliminate it .

Dear Liu,
this message most probably hints that your system is
too big to fit into memory.
Note that you don't need to use the same the same multiplication
numbers along both directions. Increase just one in order
to scan the dispersion along one direction.

Did you check that your Gamma phonons are OK,
in order to be sure what MeshCutoff you need,
and not taking it excessively high.

Best regards

Andrei

>
> sincerely
> Liu
>
> 
> NumberOfSpecies 1
> %block chemicalspecieslabel
> 1 6 C
> %endblock chemicalspecieslabel
> PAO.BasisSize DZP
> MeshCutoff   500.0 Ry/ 300Ry / 200 Ry / 100Ry
> DM.MixingWeight  0.001 / 0.01/ 0.0001
> DM.NumberPulay3
> DM.Tolerance 1.0d-4
> %include FC.fdf
>  NumberOfAtoms4
> LatticeConstant  2.45951215  Ang
> %block LatticeVectors
> 0.8660254   0.5  0.0
> 0.8660254  -0.5  0.0
> 0.0 0.0  8.0
> %endblock LatticeVectors
> AtomicCoordinatesFormat   NotScaledCartesianAng
> %block AtomicCoordinatesAndAtomicSpecies
>   -0.070.17   -0.127708   112.01
>1.4199870.12   -0.127709   112.01
>   -0.21   -0.103.725312   112.01
>1.419990   -0.053.725319   112.01
> %endblock AtomicCoordinatesAndAtomicSpecies
> SuperCell_1  4
> SuperCell_2  4
> SuperCell_3  0
> AtomicDispl  0.2  Bohr
> BandLinesScale   pi/a
> %block BandLines
> 1  0.000  0.0000.000  \Gamma
> 173  1.1547005  0.0000.000   M
> 100  1.1547005  0.7  0.000   K
> 200  0.000  0.0000.000  \Gamma
> %endblock BandLines
>
> _
> My outfile also attached. Thanks for your help again!
>
>
>
>
>
>
> At 2011-10-10 04:38:35,apost...@uni-osnabrueck.de wrote:
>>> Dear siesta users,
>>>As we konw that there should be no negative value in phonon
>>> spectrum, but when I tried to draw the phonon spectrum of
>>> bilayer-graphene, there are always some negative intervals. I'am
>>> very confused about this.
>>
>>Dear Liu,
>>
>>such kind of "noise" in the vicinity of Gamma appears
>>when the force constants are not enough converged in the real space
>>(with distance to neighbors).
>>You use the supercell (3 3 0); will (4 4 0) be affordable?
>>
>>However! - the supercell size will affect just the vicinity of Gamma,
>>but NOT Gamma itself - which ought to be good even in (0 0 0)
>>supercell. I cannot judge from your plot whether you have
>>an imaginary frequency also at Gamma, or it is just traced so
>>through many points. If you have a problem at Gamma,
>>then you should fix it first, staying with a (0 0 0) supercell.
>>It must be due to insufficient MeshCutoff
>>(170 Ry you use is not that high).
>>
>>Best regards
>>
>>Andrei Postnikov
>>
>>
>>
>>
>>>My program attached.
>>>
>>> Best regards
>>> Liu
>>>  _-
>>>
>>> NumberOfSpecies 1
>>>
>>>
>>>
>>> %block chemicalspecieslabel
>>>
>>> 1 6 C
>>>
>>> %endblock chemicalspecieslabel
>>>
>>>
>>>
>>> PAO.BasisSize DZP
>>>
>>>
>>>
>>> MeshCutoff 170.0 Ry
>>>
>>> DM.MixingWeight 0.001 /0.01/0.1
>>>
>>> DM.NumberPulay   3
>>>
>>> DM.Tolerance1.0d-5
>>>
>>> %include FC.fdf
>>>
>>> NumberOfAtoms   4
>>>
>>> LatticeConstant 2.45951215 Ang
>>>
>>>
>>>
>>> %block LatticeVectors
>>>
>>> 0.8660254  0.5 0.0
>>>
>>> 0.8660254 -0.5 0.0
>>>
>>> 0.00.0 8.0
>>>
>>> %endblock LatticeVectors
>>>
>>>
>>>
>>> AtomicCoordinatesFormat  NotScaledCartesianAng
>>>
>>>
>>>
>>> %block AtomicCoordinatesAndAtomicSpecies
>>>
>>>  -0.07   0.17  -0.127708  1   12.01
>>>
>>>   1.419987   0.12  -0.127709  1   12.01
>>>
>>>  -0.21  -0.10   3.725312  1   12.01
>>>
>>>   1.419990  -0.05   3.725319  1   12.01
>>>
>>> %endblock AtomicCoordinatesAndAtomicSpecies
>>>
>>>
>>>
>>> SuperCell_1 3
>>>
>>> SuperCell_2 3
>>>
>>> SuperCell_3 0
>>>
>>>
>>>
>>> AtomicDispl 0.2 Bohr
>>>
>>>
>>>
>>> BandLinesScale  pi/a
>>>
>>> %block BandLines
>>>
>>> 1 0.000 0.000   0.000 \Gamma
>>>
>>> 173 1.1547005 0.000   0.000  M
>>>
>>> 100 1.1547005 0.7 0.000  K
>>>
>>> 200 0.000 0.000   0.000 \Gamma
>>>
>>> %endblock BandLines
>>>
>>>
>>
>

Re: [SIESTA-L] Negative value in phonon spectrum

[apostnik]

> Dear siesta users,
>As we konw that there should be no negative value in phonon
> spectrum, but when I tried to draw the phonon spectrum of
> bilayer-graphene, there are always some negative intervals. I'am
> very confused about this.

Dear Liu,

such kind of "noise" in the vicinity of Gamma appears
when the force constants are not enough converged in the real space
(with distance to neighbors).
You use the supercell (3 3 0); will (4 4 0) be affordable?

However! - the supercell size will affect just the vicinity of Gamma,
but NOT Gamma itself - which ought to be good even in (0 0 0)
supercell. I cannot judge from your plot whether you have
an imaginary frequency also at Gamma, or it is just traced so
through many points. If you have a problem at Gamma,
then you should fix it first, staying with a (0 0 0) supercell.
It must be due to insufficient MeshCutoff
(170 Ry you use is not that high).

Best regards

Andrei Postnikov




>My program attached.
>
> Best regards
> Liu
>  _-
>
> NumberOfSpecies 1
>
>
>
> %block chemicalspecieslabel
>
> 1 6 C
>
> %endblock chemicalspecieslabel
>
>
>
> PAO.BasisSize DZP
>
>
>
> MeshCutoff 170.0 Ry
>
> DM.MixingWeight 0.001 /0.01/0.1
>
> DM.NumberPulay   3
>
> DM.Tolerance1.0d-5
>
> %include FC.fdf
>
> NumberOfAtoms   4
>
> LatticeConstant 2.45951215 Ang
>
>
>
> %block LatticeVectors
>
> 0.8660254  0.5 0.0
>
> 0.8660254 -0.5 0.0
>
> 0.00.0 8.0
>
> %endblock LatticeVectors
>
>
>
> AtomicCoordinatesFormat  NotScaledCartesianAng
>
>
>
> %block AtomicCoordinatesAndAtomicSpecies
>
>  -0.07   0.17  -0.127708  1   12.01
>
>   1.419987   0.12  -0.127709  1   12.01
>
>  -0.21  -0.10   3.725312  1   12.01
>
>   1.419990  -0.05   3.725319  1   12.01
>
> %endblock AtomicCoordinatesAndAtomicSpecies
>
>
>
> SuperCell_1 3
>
> SuperCell_2 3
>
> SuperCell_3 0
>
>
>
> AtomicDispl 0.2 Bohr
>
>
>
> BandLinesScale  pi/a
>
> %block BandLines
>
> 1 0.000 0.000   0.000 \Gamma
>
> 173 1.1547005 0.000   0.000  M
>
> 100 1.1547005 0.7 0.000  K
>
> 200 0.000 0.000   0.000 \Gamma
>
> %endblock BandLines
>
>

Re: [SIESTA-L] no PBC

[apostnik]

Sorry, I tend to disagree with Ney Henrique
that the things are that bad; please correct me if I am wrong:
If basis function of your fragments do not overlap
(i.e., the "molecule" regime is recognized) - then,
I thought, the calculation is as good as it would be
for isolated molecule in the infinite cell.
I think, notably an issue of a charged molecule
do not pose a problem and do not produce any infinite replicas,
because as the "molecule" regime is recognized, the
Madelung terms are switched off. (Note that Siesta allows
charged systems for cubic cells only, for which she knows
how to calculate Madelung sums).
I am not so sure about whether there are differences
in less sensitive issues, like construction of the average potential,
higher-order multipoles, or treatment of empty space
(no basis) on a grid... Anything else?

Best regards

Andrei Postnikov

> Hi all
>
> @ Gregorio: Yes, you are right! But strictly speaking that is not the same
> as running an isolated system. You will have infinite replicas of your
> isolated molecule anyway ... and it is problematic if you have charged
> species (there will be a rigid potential shift due to the periodic
> replication of the charges) although the Markov-Paine correction is
> implemented.
>
> @ Robert : Take a closer look to the Manual ... all this stuff is
> addressed
> there!
>
> Cheers
>
> NH
>
> 2011/10/6 Gregorio García Moreno 
>
>> **
>> Yes
>> Using supercell aproximation. i.e. Put you cluster or your aisolated
>> molecule, within of a very big cells. Thus, there are not interactions
>> between molecules of different cells,  like in a isolate molecule
>> approximation.
>>
>> El 10/6/2011 2:22 PM, robert.sed...@marge.uochb.cas.cz escribió:
>>
>> Dear Users,
>>
>> I have one probably very basic question. Is it possible to run
>> calculation
>> (for some cluster system) without periodic boundary conditions in Siesta
>> ?
>>
>> Thanks
>>
>> Robert Sedlak
>>
>>
>>
>>
>>
>
>
> --
> Dr. Ney Henrique Moreira
> Bremer Center for Computational Materials Science
> Am Fallturm 1
> 28359 Bremen
> Deutschland
> mobile: 0049176-20485882
>

Re: [SIESTA-L] Slab convergence problems

[apostnik]

Dear Carlo,
you may believe that your Mixing Weight is 0.0010
but in reality it is 0.25 (default)
because you misspelled
DM.MixingWeight
(an "s" at the end was too much; - always check the output values,
not the input ones!)

As a result, your calculation goes aghast from the 2d electronic
iteration. This has to be fixed first,
prior to MeshCutoff, k-points, wrong atom positions
(because
siesta: System type = bulk
means that you calculate not what you want),
or whatever else.

Best regards

Andrei



> Dear Andrei, Slimane, and Marcos,
>
> thank you for your suggestion.
>
>  @Andrei: you are right. I am aware of what you wrote, actually this was
> just a problem of my last run, when I forgot to delete all the output
> files
> of the previous calculation. Unfortunately, this is not the source of the
> problem.
>
> @Slimane and Marcos: I had already tried some of the hints you told me,
> but
> without success. Anyway, I stressed the values of these parameters and
> attached you can see what happens in the output file. In that calculation,
> I
> set
>
> DM.MixingWeights  0.0010
> DM.Tolerance  5E-4
> ElectronicTemperature 0.04 Ry
> MeshCutoff   400. Ry
>
> As you can see, the problem persists.
>
> Thanks again,
>
> Carlo
>
>
>
> 2011/9/20 Marcos Veríssimo Alves 
>
>> In addition to Andrei's and Slimane's suggestions, I'd increase the mesh
>> cutoff. You have too few atoms to justify such a low cutoff, and a finer
>> mesh could help in the convergence. Keep the electronic temperature high
>> from the beginning, increasing it if necessary until you reach
>> convergence
>> in the first SCF cycle, but always start the cycle from scratch (i.e.,
>> don't
>> read it from previous runs if it hasn't converged).
>>
>> Marcos
>>
>>
>> On Tue, Sep 20, 2011 at 5:37 PM,  wrote:
>>
>>> > Dear SIESTA users,
>>> >
>>> > I am doing a calculation for a Cu slab of 7 layers, exposing the
>>> (111)
>>> > surface.
>>> > I relaxing the atomic coordinates, but I have problems on the
>>> convergence
>>> > of
>>> > the scf cycles and I have no idea of the source of problem.
>>> > Attached you may find the output file.
>>> > Can someone give me a hint of a possible error?
>>>
>>> Dear Carlo,
>>>
>>> your output file says that Siesta reads previous DM from file,
>>> "fixes" it,
>>> and then immediately the numbers are crazy
>>> (already the first electronic structure iteration).
>>> Forget relaxing the coordinates.
>>> Where does this previous DM come from?
>>>
>>> In case of doubt, why don't you erase (or rename) it and start
>>> calculation from scratch?
>>>
>>> Best regards
>>>
>>> Andrei Postnikov
>>>
>>> >
>>> >
>>> > Thanks in advance,
>>> >
>>> > C. Motta
>>> >
>>>
>>>
>>
>

Re: [SIESTA-L] Slab convergence problems

[apostnik]

> Dear SIESTA users,
>
> I am doing a calculation for a Cu slab of 7 layers, exposing the (111)
> surface.
> I relaxing the atomic coordinates, but I have problems on the convergence
> of
> the scf cycles and I have no idea of the source of problem.
> Attached you may find the output file.
> Can someone give me a hint of a possible error?

Dear Carlo,

your output file says that Siesta reads previous DM from file,
"fixes" it,
and then immediately the numbers are crazy
(already the first electronic structure iteration).
Forget relaxing the coordinates.
Where does this previous DM come from?

In case of doubt, why don't you erase (or rename) it and start
calculation from scratch?

Best regards

Andrei Postnikov

>
>
> Thanks in advance,
>
> C. Motta
>

Re: [SIESTA-L] Iron 001 surface DOS

[apostnik]

Dear Deniz -
sorry, but again you don't give enough information
for giving you a meaningful advise. Let's go through it again,
step by step:

1. Forget the slab for the beginning, let's discuss bulk Fe.
Is the comparison VASP / Siesta OK? All peaks in place?
(relative to respective Fermi levels)
Magnetism OK?
The same XC used in both calculations?

2. Now you take slab, and something is shifted.
Is the DOS of the inner layer still OK?
(with respect to the Fermi level of the slab)
When does the shift starts to occur?
What says the comparison of magnetic moments
on Fe atoms, layer by layer, in VASP and Siesta?
(You shouldn't expect the numbers of local DOS to be
very close, because of their different definition in both
methods, but the layer-by-layer trends should make sense).

If you pinpoint the problem to be a sufrace effect -
then I don't know; try to make basis as complete
and long-ranged as you can;
or add ghost atoms above the surface..?
Check previous works done with Siesta on metal slabs -
not necessarily iron - for more hints...

Best regards

Andrei Postnikov


> Dear Andrei
>
> thank you for reply. I am sorry I did not write enough information. I used
> the same slab geometry for both vasp and siesta calculations.
> My k-mesh is 19x19x1 (I think it is enough for test purposes).  In both
> methods, I shifted to fermi level to  zero energy to make comparison. I am
> comparing total DOS (and also local DOS at a given layer), since i will do
> transport calculations for this iron system at the end. I am not sure that
> bulk optimized basis sets can work well for slab calculations. For this
> reason, I am trying to find good basis set for iron slab not for bulk
> iron. I attached one of possible basis set:
>
> Fe-POTCAR 3# Species label, number of
> l-shells
>  n=40   2   P
>7.0  0.
>  n=41   1
>7.0
>  n=32   2
>5.6  0.
>
> best regards,
>
> On Sep 17, 2011, at 6:29 PM, apost...@uni-osnabrueck.de wrote:
>
>>> Dear Siesta users,
>>>
>>> I am currently working on iron 001 surface.  I tried several
>>> PAO.EnergyShift and PAO.SplitNorm values to find suitable basis set.
>>> To
>>> check reliability of my basis sets   I compared siesta DOS with VASP
>>> DOS
>>> calculations.
>>
>> Dear Deniz,
>> I don't think that comparing the DOS is the best criterion
>> to decide on the quality of basis set. Anyway, from what you write
>> not much is clear. What do you compare, in fact:
>> - total DOS of you slab? (and then you make sure that you use
>> the same slab construction in both calculations)? - Or,
>> - local DOS at a given atom plane? (then you should know that
>> the definition of local DOS is completely different in two methods)
>>
>>> I observed that energy levels calculated with siesta  shift
>>> to lower energies by  0.2-0.25 eV with respect to that of VASP.
>>> What is the reason of this shift ?
>>
>> Which levels, measured how? In "absolute" values, measured
>> with respect to "zero" of energy in each method? (these two energy
>> scales are absolutely unrelated) - Or, with respect to the Fermi
>> energy of each method?
>>
>> And, your k-meshes in both methods are of course sufficiently dense
>> to be sure that the energy levels you refer to are meaningful
>> and "well converged"?
>>
>> Apart from comparing surface calculations, are you happy
>> with Siesta vs. VASP comparison for bulk iron?
>> Everything is OK with magnetic moments?
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>>
>>> Is this due to my basis set definitions
>>> or other reasons ( which i dont know)?  Could anyone supply me
>>> optimized
>>> basis set for iron surface?
>>>
>>> Thanks in advance...
>>>
>>> Deniz Cakir
>>> Computational Materials Science (CMS)
>>> Faculty of Science and Technology
>>> University of Twente
>>> Enschede, The Netherlands
>>> e-mail: d.ca...@utwente.nl
>>>
>>>
>>
>
> Deniz Cakir
> Computational Materials Science (CMS)
> Faculty of Science and Technology
> University of Twente
> Enschede, The Netherlands
> e-mail: d.ca...@utwente.nl
>
>

Re: [SIESTA-L] Iron 001 surface DOS

[apostnik]

> Dear Siesta users,
>
> I am currently working on iron 001 surface.  I tried several
> PAO.EnergyShift and PAO.SplitNorm values to find suitable basis set.  To
> check reliability of my basis sets   I compared siesta DOS with VASP DOS
> calculations.

Dear Deniz,
I don't think that comparing the DOS is the best criterion
to decide on the quality of basis set. Anyway, from what you write
not much is clear. What do you compare, in fact:
- total DOS of you slab? (and then you make sure that you use
the same slab construction in both calculations)? - Or,
- local DOS at a given atom plane? (then you should know that
the definition of local DOS is completely different in two methods)

> I observed that energy levels calculated with siesta  shift
> to lower energies by  0.2-0.25 eV with respect to that of VASP.
> What is the reason of this shift ?

Which levels, measured how? In "absolute" values, measured
with respect to "zero" of energy in each method? (these two energy
scales are absolutely unrelated) - Or, with respect to the Fermi
energy of each method?

And, your k-meshes in both methods are of course sufficiently dense
to be sure that the energy levels you refer to are meaningful
and "well converged"?

Apart from comparing surface calculations, are you happy
with Siesta vs. VASP comparison for bulk iron?
Everything is OK with magnetic moments?

Best regards

Andrei Postnikov



> Is this due to my basis set definitions
> or other reasons ( which i dont know)?  Could anyone supply me optimized
> basis set for iron surface?
>
> Thanks in advance...
>
> Deniz Cakir
> Computational Materials Science (CMS)
> Faculty of Science and Technology
> University of Twente
> Enschede, The Netherlands
> e-mail: d.ca...@utwente.nl
>
>

Re: [SIESTA-L] Stress In Wurtzite Structure

[apostnik]

> Dear Andrei,
>
> Thank you for your help. As you stated I have changed the
> AtomicCoordinatesFormat   Fractional
>
> and as Ulrich stated in the previous post, I have used
> MD.NumCGsteps 100
>
> and it worked well:

Dear Onur -
yes, of course the number of steps > 0 was crucial for getting
the system relaxed...

However:

> But the system was converging and it is converging now, so I think the
> MeshCutoff is ok (It is 350. Ry)

>From just the fact that it converges you cannot deduce anything
about whether the MeshCutoff is OK. You might get a fine convergence
with very bad MeshCutoff, or vice versa.

> Actually I thought that, I had needed to optimize PAO.EnergyShift value,
> so
> I changed it from 40 meV to 320 meV.

Beware; PAO.EnergyShift is NOT a parameter to optimize!
Because it is NOT variational, and Siesta is NOT a semiempirical method.
By increasing it, you make your basis functions shorter (more confined);
this is your free choice (hopefully based on some consideration),
similar to choosing Gaussian code other time
and planewave code yet another time.

Best regards

Andrei Postnikov

>
> Best regards.
> Onur Kavcı
>
> 2011/9/12 
>
>> Dear Onur,
>> some observations:
>>
>> 1. Your
>> AtomicCoordinatesAndAtomicSpecies
>> look like fractional ones for wurtzite, yet you have
>>  AtomicCoordinatesFormat ScaledCartesian
>> in your input file.
>> I'm afraid this combination produces not what you expect.
>> It would make sense to check (in fact, always...)
>> what is the structure as understood by Siesta
>> (and put into the .XV file).
>>
>> 2. There is a parameter
>> MD.MaxStressTol
>> which can be useful in controlling the smallness of stress
>> at which the relaxation stops. Another useful parameter
>> in this relation is
>> MD.MaxForceTol
>>
>> However, I guess this is not the main source of your problem.
>> The default MD.MaxForceTol is 0.04 eV/Ang
>> yet your forces are as large as +/- 40 eV/Ang.
>> Without knowing the details of your calculation I'd guess
>> that your structure is not converging at all. Even if started
>> from a wrong structure (see p.1), it should be able somehow
>> to converge to something. If it does not, I'd guess that
>> the MeshCutoff is the usual suspect. However, the sum of forces
>> over all atoms looks OK.
>>
>> Sorry, there is not enough information to guess further.
>> I'd only add that what you say you were doing
>> > I tried different
>> > PAO.EnergyShift values varying from 40 meV to 320 meV
>> is not a recommended way to solve this (or any other) problem.
>> PAO.EnergyShift is not a variational parameter, nor is it useful
>> to tune you results to experiment, or whatever.
>> Its role is a bit different.
>> Your another action:
>> > and I used the values
>> > 350, 270, 250, 200, 150, 100 (in Ry) for MeshCutoff.
>> Now, this parameter is not a matter of wild guess,
>> but subject to systematic test.
>> Please note that not its value as such is of importance,
>> but the number of mesh divisions it produces.
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>>
>>
>> > Dear Siesta users,
>> >
>> > I am doing some calculations about gama MnS which is in wurtzite
>> > structure.
>> > But I couldn't get the stress tensor reduced:
>> >
>> > siesta: Atomic forces (eV/Ang):
>> >> siesta:  1   -0.1430170.432695  -48.089807
>> >> siesta:  22.368858   -2.282692  -45.334125
>> >> siesta:  3   -2.0749642.005715   45.985081
>> >> siesta:  4   -0.154033   -0.154388   47.421339
>> >> siesta: 
>> >> siesta:Tot   -0.0031550.001330   -0.017512
>> >>
>> >> siesta: Stress tensor (static) (eV/Ang**3):
>> >> siesta: 0.0100290.036714   -0.017599
>> >> siesta: 0.0367120.0117960.029333
>> >> siesta:-0.0175990.029333   -1.450090
>> >>
>> >
>> > I don't know the reason of this stress tensor. I tried different
>> > PAO.EnergyShift values varying from 40 meV to 320 meV and I used the
>> > values
>> > 350, 270, 250, 200, 150, 100 (in Ry) for MeshCutoff. But the stress
>> tensor
>> > didn't decrease.  If you could suggest a solution, I would be very
>> > appreciated.
>> >
>> > I have attached the input file.
>> >
>> > Best regards.
>> > Onur Kavcı
>> >
>>
>>
>

Re: [SIESTA-L] Water vibrational frequencies

[apostnik]

;> >
>> > Thanks,
>> > Prat
>> >
>> > On Fri, Sep 9, 2011 at 2:09 AM,  wrote:
>> >
>> >> Dear Prat,
>> >>
>> >> I don't quite understand your attached report on what you call
>> >> "eggbox effect". In fact what should be expected - and can be
>> >> always enforced - is that, as MeshCutoff is increased,
>> >> the magnitude of fluctuations of total energy AND sum of forces,
>> >> as functions of uniform displacement of all atoms,
>> >> gradually decreases. This is not happening in your case.
>> >> Moreover the spatial period of fluctuations should decrease
>> >> as you increase MeshCutoff, that is not obvious in your case
>> >> (maybe, the step in displacement is too large, and the trends
>> >> are not clearly seen).
>> >> You should check that the mesh REALLY changes as you increase
>> >> the MeshCutoff value, and write down the ACTUAL numbers of mesh
>> >> divisions. It looks like not much is changing in reality
>> >> as you pass from 150 to 200 to 250 Ry.
>> >> How to make such tests were extensively discussed in tutorials.
>> >> You may have a look in my concise documentation,
>> >> http://www.home.uni-osnabrueck.de/apostnik/Lectures/SIESTA-tuto.pdf
>> >>
>> >> So the eggbox test must give you the reliable value of MeshCutoff,
>> >> and then three displacement frequencies MUST be close to zero.
>> >> In fact you have three such frequencies already
>> >> (eigenvalues 4 to 6). Suppose that once you are through with
>> >> eggbox test, you still have problematic frequencies (as those 1-3
>> >> in your present case). Then you should analyze where this comes from;
>> >> you have an eigenvector for each mode, and its imaginary frequency
>> >> tells you that as you introduce a combined displacement along such
>> >> eigenvector, the total energy goes down. This means that you are not
>> >> at equilibrium, and obviously this energy lowering cannot go on
>> >> forever; at some point you'll hit a better equilibrium. Try it.
>> >> Untill you get 6 frequencies close to zero, it doesn't make sense
>> >> do discuss the accuracy of the remaining three.
>> >>
>> >> A delicate issue may be the unit cell size, guaranteeing that
>> >> you are in the "Molecule" regime (no overlap of basis functions
>> >> across the cell boundary), and dipole-dipole interactions between
>> >> molecules are switched off. However, this may affect the accuracy
>> >> of the "meaningful" frequencies but not the fact that you should
>> >> try to get 6 zero frequencies before discussing anything else.
>> >>
>> >> Best regards
>> >>
>> >> Andrei
>> >>
>> >>
>> >> > oops had attached the wrong file previously, here is the correct
>> one.
>> >> >
>> >> > On Thu, Sep 8, 2011 at 11:27 PM, Prathibha Ramaprasad
>> >> > wrote:
>> >> >
>> >> >> Andrei,
>> >> >>
>> >> >> Attached are results from the eggbox effect for water using Javier
>> >> >> Junquera's optimised basis functions. The frequencies I get are
>> still
>> >> no
>> >> >> good.
>> >> >>
>> >> >> For 100Ry:
>> >> >> eigenvalue #   1  omega=  -743.903536842466
>> >> >>   eigenvalue #   2  omega=  -362.961235854244
>> >> >>   eigenvalue #   3  omega=  -105.075173442333
>> >> >>   eigenvalue #   4  omega= -4.529471170799031E-002
>> >> >>   eigenvalue #   5  omega=  1.353292598348921E-002
>> >> >>   eigenvalue #   6  omega=  2.296621100888934E-002
>> >> >>   eigenvalue #   7  omega=   279.287709550633
>> >> >>   eigenvalue #   8  omega=   423.852399278951
>> >> >>   eigenvalue #   9  omega=   936.394906840852
>> >> >>  Zero point energy = 0.101641 eV
>> >> >>
>> >> >> Thanks for your time and help.
>> >> >>
>> >> >> Prat
>> >> >>
>> >> >> On Wed, Sep 7, 2011 at 3:36 PM, 
>> wrote:
>> >> >>
>> >> >>> Dear Prat,
>> >> >>> 

Re: [SIESTA-L] Water vibrational frequencies

[apostnik]

> Andrei,
>
> The eggbox results I sent in the previous email were obtained using the
> script eggbox_size.sh and it seems to determine dx as pi/sqrt(meshcutoff).

Dear Prat -
sorry, I don't understand this. This might be a good script,
but I never tried it. pi is apparently unitless
and meshcutoff is in Ry; what is the meaning and units of your "dx"?
My understanding of the eggbox test is the following.
The higher MeshCutoff and larger the number of mesh divisions
along each edge of unit cell
(check these numbers; everything is in the output!),
the smaller is the length of each division, let's call it "D".
If, say, lattice vector has a length 4 Ang and number of divisions is 80,
D = 0.05 Ang.
As you uniformly displace all your atoms, you'll see fluctuations of total
energy and forces roughly with the period of "D". In fact,
in order to see it, you need SEVERAL (~10) trial displacement points
ON THE LENGTH OF "D", i.e. (in the above example)
the trial step should be about 0.005 Ang.
Then, you can make an idea of the magnitude of fluctuation.
As you increase MeshCutoff, the forces fluctuate on a smaller length,
and the magnitude of their fluctuations gets gradually suppressed,
that's is the idea.
In your figures however, I don't see any fluctuations.
The only figure which relates to forces shows smooth curves plotted
with a step of 0.1 Ang; this is too big!
Your plot is like a collection of randomly selected points
from different sinusoidal functions; nothing can be judged from it.
Other your figures show total energies (absolute values of) as functions
of the same (non indicative) displacements, plotted is such scale that
no fluctuations are seen at all. What do you conclude from all this?
Why won't you read what was written - many times earlier -
on how to do this kind of tests?

Now, apparently with a (quite small) MeshCutoff 150 Ry
you've got much better "physical" frequencies.
But still, rotational frequencies (which are supposed to be zero)
are imaginary and large, hence not good.
My suggestion: you finally make a decent
eggbox test, figure out the necessary MeshCutoff value
and bring these frequencies to zero.
If it fails, you reread my previous mail (the part concerning these
frequencies).

Good luck

Andrei

> Also I figured some of the other typos in my input. However I corrected
> these and recalculated the eggbox effect for H2O using your suggestions.
> The
> atomic forces vs origin shift shows the kind of trends you mentioned
> earlier
> except for 100 Ry meshcutoff but the total energy vs shift trends differ
> very much. I have attached these results in the attached file.
>
> Thanks,
> Prat
>
> On Fri, Sep 9, 2011 at 2:09 AM,  wrote:
>
>> Dear Prat,
>>
>> I don't quite understand your attached report on what you call
>> "eggbox effect". In fact what should be expected - and can be
>> always enforced - is that, as MeshCutoff is increased,
>> the magnitude of fluctuations of total energy AND sum of forces,
>> as functions of uniform displacement of all atoms,
>> gradually decreases. This is not happening in your case.
>> Moreover the spatial period of fluctuations should decrease
>> as you increase MeshCutoff, that is not obvious in your case
>> (maybe, the step in displacement is too large, and the trends
>> are not clearly seen).
>> You should check that the mesh REALLY changes as you increase
>> the MeshCutoff value, and write down the ACTUAL numbers of mesh
>> divisions. It looks like not much is changing in reality
>> as you pass from 150 to 200 to 250 Ry.
>> How to make such tests were extensively discussed in tutorials.
>> You may have a look in my concise documentation,
>> http://www.home.uni-osnabrueck.de/apostnik/Lectures/SIESTA-tuto.pdf
>>
>> So the eggbox test must give you the reliable value of MeshCutoff,
>> and then three displacement frequencies MUST be close to zero.
>> In fact you have three such frequencies already
>> (eigenvalues 4 to 6). Suppose that once you are through with
>> eggbox test, you still have problematic frequencies (as those 1-3
>> in your present case). Then you should analyze where this comes from;
>> you have an eigenvector for each mode, and its imaginary frequency
>> tells you that as you introduce a combined displacement along such
>> eigenvector, the total energy goes down. This means that you are not
>> at equilibrium, and obviously this energy lowering cannot go on
>> forever; at some point you'll hit a better equilibrium. Try it.
>> Untill you get 6 frequencies close to zero, it doesn't make sense
>> do discuss the accuracy of the remaining three.
>>
>> A delicate issue may be the 

Re: [SIESTA-L] Stress In Wurtzite Structure

[apostnik]

Dear Onur,
some observations:

1. Your
AtomicCoordinatesAndAtomicSpecies
look like fractional ones for wurtzite, yet you have
 AtomicCoordinatesFormat ScaledCartesian
in your input file.
I'm afraid this combination produces not what you expect.
It would make sense to check (in fact, always...)
what is the structure as understood by Siesta
(and put into the .XV file).

2. There is a parameter
MD.MaxStressTol
which can be useful in controlling the smallness of stress
at which the relaxation stops. Another useful parameter
in this relation is
MD.MaxForceTol

However, I guess this is not the main source of your problem.
The default MD.MaxForceTol is 0.04 eV/Ang
yet your forces are as large as +/- 40 eV/Ang.
Without knowing the details of your calculation I'd guess
that your structure is not converging at all. Even if started
from a wrong structure (see p.1), it should be able somehow
to converge to something. If it does not, I'd guess that
the MeshCutoff is the usual suspect. However, the sum of forces
over all atoms looks OK.

Sorry, there is not enough information to guess further.
I'd only add that what you say you were doing
> I tried different
> PAO.EnergyShift values varying from 40 meV to 320 meV
is not a recommended way to solve this (or any other) problem.
PAO.EnergyShift is not a variational parameter, nor is it useful
to tune you results to experiment, or whatever.
Its role is a bit different.
Your another action:
> and I used the values
> 350, 270, 250, 200, 150, 100 (in Ry) for MeshCutoff.
Now, this parameter is not a matter of wild guess,
but subject to systematic test.
Please note that not its value as such is of importance,
but the number of mesh divisions it produces.

Best regards

Andrei Postnikov




> Dear Siesta users,
>
> I am doing some calculations about gama MnS which is in wurtzite
> structure.
> But I couldn't get the stress tensor reduced:
>
> siesta: Atomic forces (eV/Ang):
>> siesta:  1   -0.1430170.432695  -48.089807
>> siesta:  22.368858   -2.282692  -45.334125
>> siesta:  3   -2.0749642.005715   45.985081
>> siesta:  4   -0.154033   -0.154388   47.421339
>> siesta: 
>> siesta:Tot   -0.0031550.001330   -0.017512
>>
>> siesta: Stress tensor (static) (eV/Ang**3):
>> siesta: 0.0100290.036714   -0.017599
>> siesta: 0.0367120.0117960.029333
>> siesta:-0.0175990.029333   -1.450090
>>
>
> I don't know the reason of this stress tensor. I tried different
> PAO.EnergyShift values varying from 40 meV to 320 meV and I used the
> values
> 350, 270, 250, 200, 150, 100 (in Ry) for MeshCutoff. But the stress tensor
> didn't decrease.  If you could suggest a solution, I would be very
> appreciated.
>
> I have attached the input file.
>
> Best regards.
> Onur Kavcı
>

Re: [SIESTA-L] Has anyone calculated the energy difference between diamond and graphene?

[apostnik]

> BSSE might still be important. Even though the atomic bases are the
> same, their overlap is very different. Diamond is much more compact, so
> graphene atoms would have much less complete basis because there are no
> atoms out of plane that could provide extra orbitals.
>

Dear Alexander,
you observation is valid, but - how do you suggest to implement
BSSE in practical terms for this task?
Add extra fictitious atoms between graphene planes?
(Where? How many?)
Construct a "superlattice" of which diamond and graphite will be
two different sublattices??

Best regards

Andrei

Re: [SIESTA-L] Has anyone calculated the energy difference between diamond and graphene?

[apostnik]

Dear Zhen, Ye-fei:
a new preprint addresses this problem (in a much broader context)
http://arxiv.org/abs/1109.1158
that you can use as a reference.
Their method is (at least) free from ambiguities of basis choice.

If Siesta results are too far you may wish to look more attentively
in the quality of basis set. It is important, because the basis
must be good to be able to describe two different hybridizations.
In my expectation, the LDA / GGA might give a noticeable difference
(I did not check; certainly it has been studied before).
Moreover you should check the convergence of usual cutoffs
(mesh; k-points), because you compare two different structures.
van der Waals interaction, indeed, might be
important (to be added in both systems, of course...)
On the contrary, I don't understand how BSSE may play a role
(you are the same atoms with the same basis in both sustems).

Best regards

Andrei

> Might bsse be the reason?
> Ye-fei
>
> �ҵ� iPhone
>
> �� 2011-9-9��7:22��zhuz...@msu.edu ��
>
>> Dear All,
>>
>> I am doing some test calculations on diamond and graphene. I got an
>> energy difference of 0.18 eV/atom.
>> Diamond is more stable. I think the difference is too large. Does anyone
>> has any exprirence and resultson
>> these calculation using SIESTA? I used LDA and CA. I also used same Mesh
>> cutoff, energy shift and kgrid
>> sampling (16by16by16 for diamond, 16by16by1 for graphene). I am really
>> puzzled about this. Thank you
>> very much.
>>
>> ---Zhen
>
>

Re: [SIESTA-L] Water vibrational frequencies

[apostnik]

Dear Prat,

I don't quite understand your attached report on what you call
"eggbox effect". In fact what should be expected - and can be
always enforced - is that, as MeshCutoff is increased,
the magnitude of fluctuations of total energy AND sum of forces,
as functions of uniform displacement of all atoms,
gradually decreases. This is not happening in your case.
Moreover the spatial period of fluctuations should decrease
as you increase MeshCutoff, that is not obvious in your case
(maybe, the step in displacement is too large, and the trends
are not clearly seen).
You should check that the mesh REALLY changes as you increase
the MeshCutoff value, and write down the ACTUAL numbers of mesh
divisions. It looks like not much is changing in reality
as you pass from 150 to 200 to 250 Ry.
How to make such tests were extensively discussed in tutorials.
You may have a look in my concise documentation,
http://www.home.uni-osnabrueck.de/apostnik/Lectures/SIESTA-tuto.pdf

So the eggbox test must give you the reliable value of MeshCutoff,
and then three displacement frequencies MUST be close to zero.
In fact you have three such frequencies already
(eigenvalues 4 to 6). Suppose that once you are through with
eggbox test, you still have problematic frequencies (as those 1-3
in your present case). Then you should analyze where this comes from;
you have an eigenvector for each mode, and its imaginary frequency
tells you that as you introduce a combined displacement along such
eigenvector, the total energy goes down. This means that you are not
at equilibrium, and obviously this energy lowering cannot go on
forever; at some point you'll hit a better equilibrium. Try it.
Untill you get 6 frequencies close to zero, it doesn't make sense
do discuss the accuracy of the remaining three.

A delicate issue may be the unit cell size, guaranteeing that
you are in the "Molecule" regime (no overlap of basis functions
across the cell boundary), and dipole-dipole interactions between
molecules are switched off. However, this may affect the accuracy
of the "meaningful" frequencies but not the fact that you should
try to get 6 zero frequencies before discussing anything else.

Best regards

Andrei


> oops had attached the wrong file previously, here is the correct one.
>
> On Thu, Sep 8, 2011 at 11:27 PM, Prathibha Ramaprasad
> wrote:
>
>> Andrei,
>>
>> Attached are results from the eggbox effect for water using Javier
>> Junquera's optimised basis functions. The frequencies I get are still no
>> good.
>>
>> For 100Ry:
>> eigenvalue #   1  omega=  -743.903536842466
>>   eigenvalue #   2  omega=  -362.961235854244
>>   eigenvalue #   3  omega=  -105.075173442333
>>   eigenvalue #   4  omega= -4.529471170799031E-002
>>   eigenvalue #   5  omega=  1.353292598348921E-002
>>   eigenvalue #   6  omega=  2.296621100888934E-002
>>   eigenvalue #   7  omega=   279.287709550633
>>   eigenvalue #   8  omega=   423.852399278951
>>   eigenvalue #   9  omega=   936.394906840852
>>  Zero point energy = 0.101641 eV
>>
>> Thanks for your time and help.
>>
>> Prat
>>
>> On Wed, Sep 7, 2011 at 3:36 PM,  wrote:
>>
>>> Dear Prat,
>>> judging by large imaginary frequencies, your setup is no good.
>>> It's not just a bad agreement with experiment.
>>> Apparently you take all default settings for the calculation,
>>> that is obviously not sufficient here (neither MeshCutoff nor basis).
>>> You should gradually increase MeshCutoff and make the usual eggbox
>>> test.
>>> Moreover, there are basis functions tuned by Javier Junquera
>>> for H2O on the Siesta site (under Pseudos and Bases).
>>> Why don't you try them for vibrations. Then we'll see further.
>>> Best regards
>>>
>>> Andrei Postnikov
>>>
>>>
>>> > Hi Andrei,
>>> >
>>> > These are the frequencies I get:
>>> >
>>> > eigenvalue #   1  omega=  -445.13301743729897
>>> >   eigenvalue #   2  omega=  -110.38569520939947
>>> >   eigenvalue #   3  omega=  -46.712013030918499
>>> >   eigenvalue #   4  omega= -1.12145510233318604E-002
>>> >   eigenvalue #   5  omega=  1.10614378824742104E-002
>>> >   eigenvalue #   6  omega=  1.36483674814901548E-002
>>> >   eigenvalue #   7  omega=   67.336297669281663
>>> >   eigenvalue #   8  omega=   115.87202534402807
>>> >   eigenvalue #   9  omega=   317.70618567755037
>>> >  Zero point energy =   

Re: [SIESTA-L] Water vibrational frequencies

[apostnik]

Dear Prat,
judging by large imaginary frequencies, your setup is no good.
It's not just a bad agreement with experiment.
Apparently you take all default settings for the calculation,
that is obviously not sufficient here (neither MeshCutoff nor basis).
You should gradually increase MeshCutoff and make the usual eggbox test.
Moreover, there are basis functions tuned by Javier Junquera
for H2O on the Siesta site (under Pseudos and Bases).
Why don't you try them for vibrations. Then we'll see further.
Best regards

Andrei Postnikov


> Hi Andrei,
>
> These are the frequencies I get:
>
> eigenvalue #   1  omega=  -445.13301743729897
>   eigenvalue #   2  omega=  -110.38569520939947
>   eigenvalue #   3  omega=  -46.712013030918499
>   eigenvalue #   4  omega= -1.12145510233318604E-002
>   eigenvalue #   5  omega=  1.10614378824742104E-002
>   eigenvalue #   6  omega=  1.36483674814901548E-002
>   eigenvalue #   7  omega=   67.336297669281663
>   eigenvalue #   8  omega=   115.87202534402807
>   eigenvalue #   9  omega=   317.70618567755037
>  Zero point energy = 0.031055 eV
>
> Also attached here are the input files.
>
> Thank you for your help.
>
> Prat
>
> On Tue, Sep 6, 2011 at 11:03 PM,  wrote:
>
>> Dear Prat:
>> what are your (nine) frequencies, in fact?
>> How close to zero are six of them?
>> 100 Ry cutoff a priori seems too light ...
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>> > Hi Siesta users,
>> >
>> > When I ran a vibra tool for water, I am getting ridiculously low
>> values.
>> > From what I see in the code, it seems to be converting the frequencies
>> to
>> > cm^-1 before it outputs to the file, but it looks more like the
>> numbers
>> > are
>> > in THz. I am using GGA, 100Ry mesh cutoff and 1 molecule of water in a
>> > 20X20X20 box.
>> >
>> > Has anyone here able to obtain anything close to expt values? Any
>> > pointers?
>> >
>> >
>> > Thank you,
>> >
>> > Prat
>> >
>>
>>
>

Re: [SIESTA-L] Water vibrational frequencies

[apostnik]

Dear Prat:
what are your (nine) frequencies, in fact?
How close to zero are six of them?
100 Ry cutoff a priori seems too light ...

Best regards

Andrei Postnikov


> Hi Siesta users,
>
> When I ran a vibra tool for water, I am getting ridiculously low values.
> From what I see in the code, it seems to be converting the frequencies to
> cm^-1 before it outputs to the file, but it looks more like the numbers
> are
> in THz. I am using GGA, 100Ry mesh cutoff and 1 molecule of water in a
> 20X20X20 box.
>
> Has anyone here able to obtain anything close to expt values? Any
> pointers?
>
>
> Thank you,
>
> Prat
>

Re: [SIESTA-L] Problem while plotting charge density surface

[apostnik]

> Hi
>
> I have attached an image which i had generated using the rho2xsf script
> with
> white being charge density at an isovalue of 0.04 and red being charge
> density at an isovalue of 0.08. Does this look correct?? or is there any
> mistake in it as there are abrupt discontinuities in charge densities
>

Hi
What you mean by "correct"? When you produce a figure,
you have hopefully some idea of what you see and want to emphasize
and to discuss in it.
For a piece of art, it looks fine. For scientific use,
your charge density is generated with so sparse grid
(hence the "abrupt discontinuities" you notice)
that it makes no sense to discuss it.

I think we were already through it earlier.

Best regards

Andrei Postnikov

Re: [SIESTA-L] Problem while plotting charge density surface

[apostnik]

> hi
>
> i was giving a 1 x 1 x 10 grid
> now i changed it to 10 x 10 x 20 grid
>
> and generated the xsf file

Now it's technically correct, but your grid too sparse
to see anything meaningful.

> can you please how to view the charge density surface using the Data Grid
> option in xcrysden

Select in the top menu: Tools -> Data Grid
Select the data block (there is only one in your case, already selected),
set Isovalue (in the window that opens),
press 'Submit'.

> what value should i give as isosurface value??

You see max. and min. values above the isovalue window.
Obviously you should give
a value in between, to see anything.
The rest depends on what you want to emphasize
in your figure (trial and error, please).

Best regards

Andrei Postikov

Re: [SIESTA-L] Problem while plotting charge density surface

[apostnik]

Dear Deepak,
thank you for your interest in using my scripts.
Your problems are the following:

1. When generating your RHO on a grid
(as prompted
"Enter number of grid points along three vectors"
by rho2xsf),
you apparently passed three numbers (0, 0, 10). This generated
a one-dim. grid, see the end of your "scat.XSF" file:

BEGIN_BLOCK_DATAGRID_1D
 DATA_from:scat.RHO
 BEGIN_DATAGRID_1D_RHO:spin_1
10
  0.000E+00  0.000E+00  0.000E+00
  0.000E+00  0.000E+00  1.7050300E+01
  0.0E+00  0.0E+00  0.0E+00  0.0E+00  0.0E+00 
0.0E+00
  0.0E+00  0.0E+00  0.0E+00  0.0E+00
 END_DATAGRID_1D
END_BLOCK_DATAGRID_1D

The format of this block is analogous to that of 2D and 3D blocks
in XCrySDen, but in fact XCrySDen does not support plotting 1-dim. grids,
see
http://www.xcrysden.org/doc/XSF.html#__toc__11
Which is quite understandable, because in order to plot a Y(X) function
you don't need any XCrySDen.
If this is indeed  your intention, you can use another my script,
http://www.home.uni-osnabrueck.de/apostnik/Software/rho2xy.f
created for exactly this purpose. You can then plot its result
with gnuplot, grace, or whatever.
But still, you should make a reasonable choice of your grid.
For the moment, you defined a line which runs from
(0, 0, 0) to (0, 0, 17.0503)
whereas your system contains the atoms in a plane with X=5.
So you scan parallel to your graphene at 5 Ang distance;
no wonder the result is zero everywhere.

2. What you attached under the name 'scat.LDOS' does not look like
a .LDOS file to me (which must be binary, I think) but rather
like a energy-resolved density of states ( .PDOS file?)
Again, you can visualize it with your favorite Y(X) viewer;
XCrySDen won't help you, and rho2xsf is not supposed to transform
this file.

Best regards

Andrei Postnikov


> Hi
>
> i downloaded utils from Andrei's homepage that to plot charge densities
> and
> tried to plot it using rho2xsf file. when it asks for add property to
> grid,
> if i add LDOS file it shows
>
> "At line 148 of file rho2xsf (unit=12, file = "scat.LDOS")
> Fortran runtime error; End of file"
>
> and exits
>
> however if i add the rho file i t works fine and creates the XSF file
>
> But when i open it in my xcrysden itshows error while parsing the file
>
> i have attached my ldos and XSF file with this mail
>
> --
> regards
> deepak
> 3rd Year Under Graduate
> Electronics Design and Manufacturing
> IIIT (D&M),Kancheepuram
> IIT Madras Campus,Chennai
>

Re: [SIESTA-L] experimental Dit values and the values from DFT

[apostnik]

Dear M Alaaii,

I'm afraid your explanation of what you are doing is not clear enough
for the collegues to help you; in particular you mention LDOS;
projected wavefunction; moreover, a confusion might arise from
"Siesta" definition of local density of states (i.e., space-resolved
property) from the "conventional" one (energy-resolved property
on specific atoms, i.e. Siesta's PDOS).

Having not much experience with interaces etc., I can only make
a guess of what might be problematic.
If you sum up all contributions (PDOS?)
from all interface atoms, it can indeed include much more than purely
interfacial states. First you have to identify the genuinely
interface states, both on the energy scale (looking at PDOS
and comparing it to that of pure Si and SiO2) as well as in space
(looking at LDOS calculated in the selected energy window).
This might be tricky because of overlap of these state with
others. However, an attempt to isolate these states
in both energy and spacial domains might give you the best you may
probably achieve. Then, you interpolate the LDOS constructed
in carefully selected energy interval onto the carefully selected
spatial box that confines possibly clearly the interface states
(using my rho2xsf script or other).
Then, as you get the desired property of a selected 3-dim. grid,
you can do with it whatever you want - sum up over Z and get it
(X,Y)-resolved; sum it up over all box and attribute to surface unit, etc.
Of course I don't know whether this will give a good agreement with expt.,
because you should look into what exactly the experiment measures.
Its reported number is possibly also obtained in some indirect way...

Best regards

Andrei Postnikov


> Dear All,
> This is my second email and I highly appreciate your help in this regard.
> I am working on Si/SiO2 system, at the moment.
>
>   1. I was wondering whether the sum of all LDOS calculated from projected
>   wavefunction on the Si at the interface can give me the density of
> interface
>   states. The question is whether the sum of LDOS of projected
> wavefunction
> on
>   the inter-facial atoms is equivalent to planar local density of state at
> the
>   interface.
>   2. If positive, the sum of LDOS of all projected wavefunction on all
>   atoms at the interface divided by the surface of the interface is by far
>   higher than the Dit reported from experiment. Actually, I am thinking,
> no
>   matter how big the structure is, the Dit is still very high. I was
> wondering
>   whether I miss any point to consider for calculation.
>
> I really really appreciate your help.
> Have a great summer.
> M Alaaii
>

Re: [SIESTA-L] got a question about meshcutoff vs.convergence

[apostnik]

>
> To control the predefined "meaningful energy DIFFERENCES", which
> parameter(s)/flag(s) you refer to in the fdf file?
>
> Regards

This is not in the flags; this is human's part of job
("to know what you are doing").

To give you an example:
when comparing energies of two different phases
we'd probably deal with another scale of
"meaningful energy DIFFERENCE"
than when comparing two different orientations of a given spin
(although who knows...)

Best regards

Andrei Postnikov



> On 08/04/2011 03:59 AM, apost...@uni-osnabrueck.de wrote:
>>> On 08/04/2011 02:17 AM, apost...@uni-osnabrueck.de wrote:
 2. You can eventually extract total energies, make a plot and
 look at it, but attention! It is IMPOSSIBLE to decide whether the
 energy
 is "converged" or not without specifying the relevant order of
 magnitude
 of DIFFERENCES you are interested in to study, like e.g. the effect of
 atom displacement, changing magnetic configuration, imposing pressure,
 etc.
>>>
>>> For total energies, this can be controlled by
>>> DM.Require.Energy.Convergence and DM.Require.Energy.Convergence.  But
>>> in
>>> fact, in default case, the siesta don't use this convergence criteria.
>>> Regards
>>> --
>>> Hongyi Zhao
>>
>> This seems to be a misunderstanding. The convergence entries you refer
>> to,
>> along with some others, determine when the program converges
>> after a certain number of iterations and stops
>> FOR THE GIVEN PARAMETERS as MeshCutoff, k-mesh density, basis size, etc.
>> This does not give you any hint as to whether these parameters
>> are "good". Whatever was referred to before (and in the mailing list,
>> and in the tutorials) in the present context as convergence with respect
>> to these calculation parameters is different.
>> What was meant is the following:
>> you take some MeshCutoff, converge the calculation,
>> write down the value. You change MeshCutoff, you converge the
>> calculation,
>> write down the value. Then you have a look at a sequence of converged
>> values you wrote down and try to figure out which value of MeshCutoff
>> is sufficient, or "converged", in a sense that its further changing
>> does not affect the result. The problem hereby: you cannot decide on
>> this
>> unless you know IN ADVANCE which accuracy you in fact need, i.e.,
>> meaningful energy DIFFERENCES for your problem.
>
> To control the predefined "meaningful energy DIFFERENCES", which
> parameter(s)/flag(s) you refer to in the fdf file?
>
> Regards
>> This was the meaning of my comment.
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>>
>>
>
>
> --
> Hongyi Zhao 
> Institute of Semiconductors, Chinese Academy of Sciences
> GnuPG DSA: 0xD108493
>
>

Re: [SIESTA-L] Energy bands and basis functions

[apostnik]

> Dear SIESTA users,
>
> Is there a straightforward way to find out which band in energy
> dispersion, E(k), relates to which basis function?

An energy band at a given k-point is an EIGENVALUE
whose corresponding EIGENVECTOR exactly shows the weight of each
basis functions in it. Basis functions are numbered throughout
according to the list of atoms in the system, as the basis functions
(according to atom's type) are merged into the total list,
atom by atom.
I think if you demand printing out the wavefunctions you'll get
this information more or less straightforwardly.
In general, EACH band includes some contribution from ALL basis functions.

Best regards

Andrei Postnikov

Re: [SIESTA-L] got a question about meshcutoff vs.convergence

[apostnik]

> On 08/04/2011 02:17 AM, apost...@uni-osnabrueck.de wrote:
>> 2. You can eventually extract total energies, make a plot and
>> look at it, but attention! It is IMPOSSIBLE to decide whether the energy
>> is "converged" or not without specifying the relevant order of magnitude
>> of DIFFERENCES you are interested in to study, like e.g. the effect of
>> atom displacement, changing magnetic configuration, imposing pressure,
>> etc.
>
> For total energies, this can be controlled by
> DM.Require.Energy.Convergence and DM.Require.Energy.Convergence.  But in
> fact, in default case, the siesta don't use this convergence criteria.
> Regards
> --
> Hongyi Zhao 

This seems to be a misunderstanding. The convergence entries you refer to,
along with some others, determine when the program converges
after a certain number of iterations and stops
FOR THE GIVEN PARAMETERS as MeshCutoff, k-mesh density, basis size, etc.
This does not give you any hint as to whether these parameters
are "good". Whatever was referred to before (and in the mailing list,
and in the tutorials) in the present context as convergence with respect
to these calculation parameters is different.
What was meant is the following:
you take some MeshCutoff, converge the calculation,
write down the value. You change MeshCutoff, you converge the calculation,
write down the value. Then you have a look at a sequence of converged
values you wrote down and try to figure out which value of MeshCutoff
is sufficient, or "converged", in a sense that its further changing
does not affect the result. The problem hereby: you cannot decide on this
unless you know IN ADVANCE which accuracy you in fact need, i.e.,
meaningful energy DIFFERENCES for your problem.
This was the meaning of my comment.

Best regards

Andrei Postnikov

Re: [SIESTA-L] got a question about meshcutoff vs.convergence

[apostnik]

Dear Lily,
I don't know what was the error in your attempted extracting lines
from files; this is not Siesta but Unix stuff; most probably you messed up
with >, <, and spaces (or, lines you searched for were not in your files).
However, I'd have couple of comments on the meaning of the test
you are doing.

1. Not the input MeshCutoff matters, but the actual numbers
of division along lattice vectors accepted by the system.
Search in the output for e.g.:
| InitMesh: MESH =40 x40 x40 =   64000
| InitMesh: Mesh cutoff (required, used) =   350.000   425.493 Ry
Therefore, small changes of input MeshCutoff might not change the mesh
at all.

2. You can eventually extract total energies, make a plot and
look at it, but attention! It is IMPOSSIBLE to decide whether the energy
is "converged" or not without specifying the relevant order of magnitude
of DIFFERENCES you are interested in to study, like e.g. the effect of
atom displacement, changing magnetic configuration, imposing pressure, etc.


Best regards

Andrei


> Dear All,
>
> When i tried to look into how the meshcutoff  affects the total energy,  i
> used several meshcutoff values, if I take CH4 as an example,   25, 35, 55,
> 85, and 125,
>
> to avoid manually recording energy and meshcutoff,   I followed the H2O
> example, in which they investigate how cell size influences the energy
>
> using commands like:
>
> siestah2o.5ang.out &tail -f h2o.5ang.out
> ...  
>
> siestah2o.35ang.out &tail -f h2o.35ang.out
>
> grep " Total=" h2o.*ang.out >convergencecellsize.dat
>
>
> then you will obtain your size and energy in each case in
> convergencecellsize.dat file.
>
> When I used
> siesta ch4.25mesh.out & tail -f ch4.25mesh.out
> ...
> siesta ch4.125mesh.out &tail -f ch4.125mesh.out
>
> grep " Total=" ch4.*mesh.out >convergencemesh.dat
>
> Somehow, convergencemesh.dat  is empty in my case,  what's wrong?
>
> Many many thanks for any hint!
>
> Lily
>

Re: [SIESTA-L] a naive question about local density of states

[apostnik]

> Dear all,
>
> when we instruct to write the LDOS, we need give two energies,( in the
> manual, it says the energies of the range for LDOS integration, relative
> to
> the program's zero, i.e. the same as the eigenvalues printed by the
> program)   What does this mean?

The energy scale in question is the same as appears in DOS / PDOS

> When I look up the output file, there
> are only 3 blocks which contain energy related information,
> the first one is the  in each iteration, there are some Eharris and E_KS
> the second one is about the program energy decomposition , and the third
> one is the final energy.

None of these... These are TOTAL energies

> Where are the eigenvalues?

In the .EIG file

> what it means " the energies of the range for LDOS integration,
> relative to the program's zero,"

A number in the relevant scale which is printed out in .out file is
Ef(eV). If it is say -4.5 eV and you are interested inintegrating LDOS
through 2 eV of the highest occupied states you define the energy window
[ -6.5 : -4.5 ]

Best regards

Andrei Postnikov

Re: [SIESTA-L] About the imaginary part of wavefunction at gamma point.

[apostnik]

> On 07/16/2011 10:58 AM, Hongyi Zhao wrote:
>> On 07/15/2011 09:13 PM, apost...@uni-osnabrueck.de wrote:
 Hi all,

 When I do the calculation on graphene, I let it write out the
 wavefunction and find that:

 For Gamma point, the value of Im(psi) always equals to zero. But I
 cann't understand why it should like this? Any hints?
>>>
>>> Dear Hongyi
>>> I think, it SHOULD not but CAN be chosen that way,
>>> and that's what the program makes by default,
>>> whereas for general k value it CANNOT,
>>> due to Bloch's periodic form of wave function:
>>> [Bloch amplitude, periodic] * exp(ikr)
>>
>> Do you mean we can select special K point as the original point of the
>> reciprocal lattice to meet the above requirement. If so, does this point
>> exist and only exist for a specific shystem?
>
> On second thought, according to the periodic form of wave function,
> there should be infinite number of such k points which meet the above
> requirement.
>

Dear Hongyi,
I cannot follow your argumentation - sorry - and cannot understand
what in your opinion I mean... My point is simply that
exp(ikr) with k=0 is real
whereas at a general k WITHIN THE BRILLOUIN ZONE, i.e. NOT pi*N,
it is complex.
(It seems however that the wave function may be chosen real
at some symmetric points at the Brillouin zone boundary, as well,
but this was not your original question...
And the number of such points is not infinite, anyway).

Have a nice day

Andrei Postnikov

Re: [SIESTA-L] question about post-process output file .PDOS

[apostnik]

> Dear all,
>
> How to give those parameters which show up as the head of .PDOS file?
>
>  pdos>
> 1
>6
> 
>
> It seems that we need modify them, otherwise, it gives us the error
> message
> like:
>
> At line 53 of file eig2dos.f (unit = 5, file = 'stdin')
> Fortran runtime error: Bad real number in item 1 of list input
>
> Any hint would be much appreciated!
>
> Lily
>

Dear Lily -
if I remember correctly,
eig2dos utility is intended to read .EIG
and has nothing to do with .PDOS
In case of doubt,
read comments at the beginning of eig2dos.f

Best regards

Andrei Postnikov

Re: [SIESTA-L] About the imaginary part of wavefunction at gamma point.

[apostnik]

> Hi all,
>
> When I do the calculation on graphene, I let it write out the
> wavefunction and find that:
>
> For Gamma point,  the value of Im(psi) always equals to zero. But I
> cann't understand why it should like this?  Any hints?

Dear Hongyi
I think, it SHOULD not but CAN be chosen that way,
and that's what the program makes by default,
whereas for general k value it CANNOT,
due to Bloch's periodic form of wave function:
[Bloch amplitude, periodic] * exp(ikr)

And certainly that's not only for graphene :-)

Best regards

Andrei Postnikov

> Regards
> --
> Hongyi Zhao 
> Institute of Semiconductors, Chinese Academy of Sciences
> GnuPG DSA: 0xD108493
>
>

Re: [SIESTA-L] pdos

[apostnik]

> Hi,
> SOrry that I didn`t know your name. I just wrote what in e-mail. and I
> thought
> may be it is your nickname or some things. any way sorry.

Hi,
nothing big to be sorry about - the anonymous
>> From: root 
mail was not from me but from another person who
does not sign by name, advertises somebody else's -
not at all anonymous - product and gives wrong advises
how to use it. I do not care much about credit for
my small script, but find such attitude somehow not correct.
I hope you can successfully use it now as you managed to compile it.

Best regards

Andrei Postnikov


>
> thanks for your complet and step by step explenation. I just compile the
> fmpdos.f, by f95, I had to typed,
>
> f95 -o name of the file after compiling  fmpdos.f
> and it compiled without any errore. now I am going to use it on my pdos to
> get
> results and by your complet explenation I am sure I can do that.
> Thanks for your wonderfull guids again,
> Good luck
>
>
>
>
> 
> From: "apost...@uni-osnabrueck.de" 
> To: siesta-l@uam.es
> Sent: Sat, July 2, 2011 1:59:58 PM
> Subject: Re: [SIESTA-L] pdos
>
>> hi dear root
>> I got fmpdos.f from internet.
>
> Dear Zahra -
> sorry for intervening even as I am not "dear root".
> The internet is big. I am glad that you find my small script
> worth using. It was intended to be simple in use;
> I am sorry that it resulted in difficulties.
>
>> but I don`t know how I have to compile it,
>> and how
>> run it and get results from my case.pdos. my compiler is f95.
>
> How about trying
> f95 fmpdos.f
> as a first guess
> (to compile) and then
> ./a.out
> (to run)? For more advanced solutions, consult your Fortran manual.
>
>> I don`t know how to insert n, m and l.
>
> for each of them, you type a number when prompted on the command line
> and press "Enter" (right on the keyboard, with big arrow on it).
>
>> dos the program ask for them when I
>> compile it.
>
> Not when you compile, but when you run the compiled code.
>
>> And I cannot understand why you wrote that l,m are always 0. in the way
>> that I
>> know the l for orbital p is equal to 1.
>
> They are not always 0, you must have misunderstood
> the explanations of dear root.
> The script prints on your screen (when started):
> ' Extract data for atom index',
> ' (enter atom NUMBER, or 0 to select all),'
> ' or for all atoms of given species',
> ' (enter its chemical LABEL):'
> to which you type, say,
> 178
> (to select the DOS on your atom Nr 178)
> or, say
> Md
> (to sum up the DOS over all Md atoms present in your system).
> As concerns the "l" value - you are right, it is indeed 1 for p orbitals.
> However,
> "n" , when prompted as
> ' Extract data for n= ... (0 for all n ):'
> , is not the "certain atom number" as dear root explains you,
> but a principal quantum number of the orbital in question.
>
> The rest, I hope, is self-evident.
>
> Good luck
>
> Andrei Postnikov
>
>
>>
>> 
>> From: root 
>> To: siesta-l@uam.es
>> Sent: Wed, June 29, 2011 11:58:11 PM
>> Subject: Re: [SIESTA-L] pdos
>>
>> Hi Zahra,
>>
>> There is a very handy script fmpdos.f.
>>
>> Follow the instruction, enter value of n, l, m, (n define certain atom
>> number,
>> l,m enter 0 for all case) define the name of output file, eg. atom 4 I
>> define it
>> as 4.outAnd you can get the DOS for any atoms you want.
>>
>>
>>
>>
>> On Tue, Jun 28, 2011 at 1:40 AM, Zahra Talebi 
>> wrote:
>>
>> Hi every body
>>>I want to draw a diagram of my case.PDOS out put after running siesta
>>> and
>>> I
>>>don`t know how to do that.
>>>
>>>If any body knows let me know step by step
>>>
>>>
>>

Re: [SIESTA-L] pdos

[apostnik]

> hi dear root
> I got fmpdos.f from internet.

Dear Zahra -
sorry for intervening even as I am not "dear root".
The internet is big. I am glad that you find my small script
worth using. It was intended to be simple in use;
I am sorry that it resulted in difficulties.

> but I don`t know how I have to compile it,
> and how
> run it and get results from my case.pdos. my compiler is f95.

How about trying
f95 fmpdos.f
as a first guess
(to compile) and then
./a.out
(to run)? For more advanced solutions, consult your Fortran manual.

> I don`t know how to insert n, m and l.

for each of them, you type a number when prompted on the command line
and press "Enter" (right on the keyboard, with big arrow on it).

> dos the program ask for them when I
> compile it.

Not when you compile, but when you run the compiled code.

> And I cannot understand why you wrote that l,m are always 0. in the way
> that I
> know the l for orbital p is equal to 1.

They are not always 0, you must have misunderstood
the explanations of dear root.
The script prints on your screen (when started):
 ' Extract data for atom index',
 ' (enter atom NUMBER, or 0 to select all),'
 ' or for all atoms of given species',
 ' (enter its chemical LABEL):'
to which you type, say,
178
(to select the DOS on your atom Nr 178)
or, say
Md
(to sum up the DOS over all Md atoms present in your system).
As concerns the "l" value - you are right, it is indeed 1 for p orbitals.
However,
"n" , when prompted as
' Extract data for n= ... (0 for all n ):'
, is not the "certain atom number" as dear root explains you,
but a principal quantum number of the orbital in question.

The rest, I hope, is self-evident.

Good luck

Andrei Postnikov


>
> 
> From: root 
> To: siesta-l@uam.es
> Sent: Wed, June 29, 2011 11:58:11 PM
> Subject: Re: [SIESTA-L] pdos
>
> Hi Zahra,
>
> There is a very handy script fmpdos.f.
>
> Follow the instruction, enter value of n, l, m, (n define certain atom
> number,
> l,m enter 0 for all case) define the name of output file, eg. atom 4 I
> define it
> as 4.outAnd you can get the DOS for any atoms you want.
>
>
>
>
> On Tue, Jun 28, 2011 at 1:40 AM, Zahra Talebi 
> wrote:
>
> Hi every body
>>I want to draw a diagram of my case.PDOS out put after running siesta and
>> I
>>don`t know how to do that.
>>
>>If any body knows let me know step by step
>>
>>
>

Re: [SIESTA-L] Liquid metals

[apostnik]

>> dear   siesta users,
>> can we use siesta to study liquid metals? If it is, how?
>>
>
> Dear Surjeet,
>
> It depends on the properties you want to study and what computational
> resources you have at your disposal. Actually we performed calculations
> of the structure of liquid rubidium, see
> http://dx.doi.org/10.1088/0953-8984/20/11/114104 for details. I know as
> well, that Dr. Postnikov et al performed calculations of liquid metals
> conductivity using Kubo-Greenwood formalism, but I don't have the
> reference right now.

True; thanks Andrey to remind about it.
Here is the reference:
F. Knider, J. Hugel and A. V. Postnikov.
Ab initio calculation of dc resistivity in liquid Al, Na and Pb.
J. Phys.: Cond. Matter 19(19), 196105 (May 2007)

Best regards

Andrei Postnikov

Re: [SIESTA-L] Question about Denchar

[apostnik]

Dear Artem,
some short comments:

1. You are right, it is a very good idea FIRST get some idea
about band structure, well converge things, AND THEN
select k-points and functions for plotting, not other way around.
2. In k-points other than Gamma, the wave functions are complex;
you should figure out how to plot them in a meaningful way;
it may seem cumbersome, though.
3. The number of occupied states in a nonmagnetic dielectric equals
(number of valence electrons)/2, i.e. exactly given by chemical
composition, prior to any calculation.
If magnetic moment is known, getting HOMO and LUMO for each spin channel
is equally straightforward. In a metal, this is of course more tricky,
as values from one k-point to another.

Best regards

Andrei Postnikov


> Dear Andrei,
>
> Thank you very much for your assistance. It seems that I overlooked the
> option that you pointed out.
>
> But, still there is a problem. WaveFuncKPoints block does allow to specify
> the k-points at which wavefunctions are calculated. This is how we can
> control this aspect at the stage of main Siesta program. But what if I
> simply do not know the whole band structure before running Siesta which
> means that I do not know the numbers of wavefunctions that I am
> potentially
> interested in (e.g. I am interested in the bands in vicinity of Fermi
> level
> but I do not know its value and numbers and values of the corresponding
> eigenstates). Also, before having the whole band structure calculated it
> is
> hard to determine in which k-points I want the wavefunctions to be
> calculated.
> For this reason, I thought that we could first let Siesta calculate the
> whole picture, then visualize the Band Structure and only after that
> decide
> which k-points and which eigenstates we are interested in. In other words,
> we could perform the control at the stage of using Denchar utility without
> reducing the outcome of calculations from Siesta. Does it make sense?
>
> Besides this issue I am not quite sure that I understand the meaning of
> option NumberOfEigenStates (*integer*). There is a restriction on that
> integer number (greater than the number of occupied states and less than
> the
> number of basis functions). How can I know the number of occupied states
> (even at T=0) before the calculation of Fermi level and "HOMO" band are
> done? I suppose that this question is related to the above mentioned.
>
> I would be very obliged if you could explain this aspect and help me with
> my
> initial problem.
>
> Best regards,
>
> Artem
>
>
> On Sat, May 28, 2011 at 3:52 AM,  wrote:
>
>>
>> >
>> > Now MY QUESTION: how to control k-point in which I want the
>> wavefunctions
>> > to
>> > be calculated, and also how to control the number of eigenstate (fist,
>> > 31st,
>> > etc), for a given k-point, for which I want the wavefunction to be
>> > calculated?  Does anybody know how to perform that controlling?
>>
>> Dear Artem,
>> I think, the description of the block
>>
>> WaveFuncKPoints
>>
>> contains the information you searches?
>> (k-points AND wavefunc. numbers to store)
>>
>> Best regards
>>
>> Andrei Postnikov
>>
>> >
>> > Many thanks for any help and your time,
>> >
>> > Artem
>> >
>>
>>
>

Re: [SIESTA-L] Question about Denchar

[apostnik]

>
> Now MY QUESTION: how to control k-point in which I want the wavefunctions
> to
> be calculated, and also how to control the number of eigenstate (fist,
> 31st,
> etc), for a given k-point, for which I want the wavefunction to be
> calculated?  Does anybody know how to perform that controlling?

Dear Artem,
I think, the description of the block

WaveFuncKPoints

contains the information you searches?
(k-points AND wavefunc. numbers to store)

Best regards

Andrei Postnikov

>
> Many thanks for any help and your time,
>
> Artem
>

Re: [SIESTA-L] is this result reliable?

[apostnik]

> Respected Prof. Andrei Postnikov and all siesta users,
>
> I have few queries..
>
>>
>> However, you should know what you are doing, because the structure
>> of wave atomic function inside the cutoff radius indeed starts to
>> differ,
>> and the error increases as you come closer.
>>
>
>> Instead of paying to much attention to wave-functions, shouldn't we
>> attend
> to norm-conservation condition (to potentials generated by ATOM code in
> siesta). Because satisfying norm-conservation (equal charge density inside
> cutoff radii for both all electron and pseudo-wave functions) implies that
> logarithmic derivatives and first energy derivatives of logarithmic
> derivatives
> corresponding to all electron and pseudo-wave functions agree at cutoff
> radii.
> Which means better scattering properties of the ion core.
>
> or may be this condition is best satisfied only by best choice of cut-off
> radii for
> each angular momentum components at particular energy or at multiple
> energy
> references.
>

Sorry Sonu,
we enter here a sensitive issues of pseudopotential construction;
I don't think I am a right person to give you an enlightening answer...


>> Moreover I think this is basically difficult to cover in Siesta,
>> using the same fixed basis, a large range of varying interatomic
>> distances
>> with the same accuracy.
>
> since the basis is fixed, but if i increase the size  of basis  say DZ -->
> DZP ---> TZP ---> TZP and diffuse functions, can i expect better behaviour
> of basis over a range of interatomic distances, say from 0.9 Ang to 2.0
> Ang
> in system under consideration ?

The short answer: yes you can. However, a priori you don't know how exactly.
This is the standard problem of fixed bases. As you, say, change volume
and your wave functions get more compressed, with the planewave basis,
or with numerical adjustable basis (LMTO, LAPW) this will be taken care of
automatically. Whereas with a fixed basis, you basis quality may be
good either here, or there, but hardly everywhere. Unless your basis
is "very" large, or you were very smart in constructing it.
But then it becomes an issue of tests and human decision
which (in)accuracy you are ready to tolerate, in your particular problem.

>> > If not, shall I reduce the cutoff radii of my input file for the
>> > pseudopotential generation?
>>
>> This is good for transferability but results in harder pseudopotential
>> and introduce other kind of problems.
>>
>> will you please elaborate, what kind of problems in addition to
> computational
> time ?

Weirdly behaving pseudopotential and pseudofunctions ->
their long-range extension in reciprocal space ->
high cutoff needed -> memory and time... Nothing more, indeed

Best regards

Andrei


>> Thank you very much.
>
> With regards,
>
> Sonu Kumar
> Phd Student
> IITD
>

Re: [SIESTA-L] siesta: WARNING: Atoms 460 641 too close: rij = 0.351881 Ang

[apostnik]

> Yeah, Usually how to come up with reasonable 'initial guess' of
> coordinates
> instead of wild guessing?
>
> For example, carbon tube, I used one ready-made program which generated
> coordinates. But it gave me the warning like " WARNING: Atoms
>>> 460 641 too close: rij = 0.351881 Ang".
>
> Many thanks
>
> Lily

Dear Lily,
try any visualization code to have a look at your generated structure
(e.g., use xv2xsf from my Sies2xsf tools, then XCrySDen)
and at the problem atoms, then you'll have a better idea what's going on.
Typical sources of problem: wrong length unit (lattice constant),
wrong translation vectors (the atoms are translated across the box
boundary where they are not supposed to be), or human typo.

Best regards

Andrei


>
> On Sun, May 15, 2011 at 12:38 AM, sonu kumar <1009uku...@gmail.com> wrote:
>
>> Hi all Siesta users,
>>
>>
>>
>> >> how to fix the question like"siesta: WARNING: Atoms
>> >> 460 641 too close: rij = 0.351881 Ang" ?
>>
>> This warning comes because you have incorrectly specified the
>> positions of atoms in your unit cell.
>>
>>  regards,
>>
>> --
>> Sonu Kumar
>> Phd Student
>> Physics Department
>> Indian Institute of Technology
>> Delhi-110016, India
>> web:-http://www.iitd.ac.in/
>>
>

Re: [SIESTA-L] is this result reliable?

[apostnik]

> Dear all,
>
> I want to do a test calculation on the local magnetic moment of iron on
> the
> interatomic distance. The cutoff radii for my pseudopotential generation
> is
> 2.0 Bohr. When the Fe-Fe distance is set to be 3.0 Bohr (this means that
> one
> atom is in the core region of the other atom), is the calculation
> reliable?

Dear Leila,

I remember to have posed such question long ago to one of Siesta fathers
and the answer was positive, in a sense that there is no constraint
technically forbidding overlap of pseudopotential cutoff radii
at adjacent atoms, like say muffin-tins in FLAPW codes.
(More specific comments from pseudopotential gurus will be mostly welcome).
However, you should know what you are doing, because the structure
of wave atomic function inside the cutoff radius indeed starts to differ,
and the error increases as you come closer.

Discussing specifically Fe, the overlap of 3d will probably not
pose big problem, because their nodal structure is correct.
For 4s, it is not, but they are delocalized anyway, and I don't think
the details of their pseudopotential count too much.
However, you should take care of 3p, including them
in the valence as you come for closer distances.

Moreover I think this is basically difficult to cover in Siesta,
using the same fixed basis, a large range of varying interatomic distances
with the same accuracy. That is, you should check against all-electron
calculation whether the error is acceptable for your purposes.

> If not, shall I reduce the cutoff radii of my input file for the
> pseudopotential generation?

This is good for transferability but results in harder pseudopotential
and introduce other kind of problems.

> Or, can I do all-electronic calculations with
> SIESTA?

No

Best regards

Andrei

>
> Thank you in advance!
>
> Best wishes!
>
> Yours sincerely,
>
> Leila
>

Re: [SIESTA-L] siesta: WARNING: Atoms 460 641 too close: rij = 0.351881 Ang

[apostnik]

> Yeah, Usually how to come up with reasonable 'initial guess' of
> coordinates
> instead of wild guessing?
>
> For example, carbon tube, I used one ready-made program which generated
> coordinates. But it gave me the warning like " WARNING: Atoms
>>> 460 641 too close: rij = 0.351881 Ang".

Dear Lily,

nobody knows where your atoms 460 and 641 are. But, using
any visualization program and looking at generated structure
with your problematic atoms will immediately give you a hint
whom to blame: the atoms, the ready-made program, or your own typo.
Typical errors are wrong length scale (lattice parameter),
or wrong translation vectors.

Best regards

Andrei



>
> Many thanks
>
>
> Lily
>
> On Sun, May 15, 2011 at 12:38 AM, sonu kumar <1009uku...@gmail.com> wrote:
>
>> Hi all Siesta users,
>>
>>
>>
>> >> how to fix the question like"siesta: WARNING: Atoms
>> >> 460 641 too close: rij = 0.351881 Ang" ?
>>
>> This warning comes because you have incorrectly specified the
>> positions of atoms in your unit cell.
>>
>>  regards,
>>
>> --
>> Sonu Kumar
>> Phd Student
>> Physics Department
>> Indian Institute of Technology
>> Delhi-110016, India
>> web:-http://www.iitd.ac.in/
>>
>

Re: [SIESTA-L] how to find bulk modulus for hexagonal cell

[apostnik]

Thank you Roberto,
that's interesting. Here is some bunch of formulae:
http://sepwww.stanford.edu/data/media/public/docs/sep125/jim1/paper_html/node10.html

Still, I don't understand the formula
P = B \delta V / V
in two aspects: 1) what is \delta? A difference? Between what and what?
2) which pressure (uniaxial? hydrostatic? some mixture?) is to monitor
as one varies the volume, keeping c/a.
About inner degrees of freedom in hcp 2-atom cell I still don't understand
either, nor about "distortion of the basal plane" as one uniformly scales
all the cell dimensions.
Whatever; this is far from the Akshu's original question.

Best regards,

Andrei

>
>  Hi Andrei and Akshu
>
>> My excuses,
>> I strongly disagree with Roberto's opinion:
>>
>
> Well, not that much ;-)
>
>>>  The closest thing to bulk modulus in hcp is obtained by varying
>>>  all the cell dimensions simultaneously and homogeneously and
>>>  monitoring the pressure while doing that. Then P = B \delta V / V.
>>
>> The definition of bulk modulus (B) known to me is
>> B = -V dP/dV.
>> As the structure is not cubic, varying cell dimensions simultaneously
>> will introduce stress (the relaxed cell won't like to have the same
>> c/a for different volumes). I think the right thing to do is
>> to vary TARGET PRESSURE, let cell parameters free and monitor the
>> VOLUME,
>> then use the formula above.
>> Or - probably better - you extract elastic constants,
>> checking back their definition in hcp, as indepent parameters.
>> Because they may behave differently; the bulk modulus is just
>> some averaged combination of them, and won't tell you much,
>> for a serious test.
>>
>
>  The formula I quoted, P = B \delta V / V,  will give you just that
>  average of elastic modulii Andrei's mentioning
> B = 1/9 {2C11+2C12+4C13+C33}
>  This, I seem to recall, is called Voigt's average. Namely one where
>  an average deformation is imposed.
>
>  Andrei's point of view is an alternative possibility that corresponds
>  to impose an average stress  (currently a pressure, Reuss average).
>
>  Clearly, the former method requires fixed cell while the latter
>  variable cell. I believe, however, that my method is a bit better
>  for numerical calculations, because it is easier (or more exact) to
>  control lattice dimensions than pressure.
>
>
 should i fix the atomic positions?
>>>
>>>  No. The hcp 2-atoms cell possesses inner degrees of freedom,
>>
>> How's that??? I always thought there is something like
>> (0 0 0) and (1/3  2/3  1/2), no internal coordinates.
>> These relative coordinates should not change
>> (they will in fact, slightly, due to lack of symmetry constraint
>> in Siesta).
>>
>
>  A distortion of the basal plane will generally  produce a relative
>  displacement between the two atoms of the cell.
>  This is because the atomic positions are not centrosymmetric.
>  I agree though, that for a uniform expansion/contraction this will
>  not happen.
>
>  Best,
>
>  Roberto
>
>

Re: [SIESTA-L] how to find bulk modulus for hexagonal cell

[apostnik]

My excuses,
I strongly disagree with Roberto's opinion:

>
>  Hi akshu,
>
>> dear siesta users
>> i need to calculate the bulk modulus in order to test the
>> pseudopotential
>> for technetium. Tc is found in hexagonal structure, so along with
>> varying
>> the lattice constant, do i need to change 'c' as well?
>
>  The closest thing to bulk modulus in hcp is obtained by varying
>  all the cell dimensions simultaneously and homogeneously and
>  monitoring the pressure while doing that. Then P = B \delta V / V.

The definition of bulk modulus (B) known to me is
B = -V dP/dV.
As the structure is not cubic, varying cell dimensions simultaneously
will introduce stress (the relaxed cell won't like to have the same
c/a for different volumes). I think the right thing to do is
to vary TARGET PRESSURE, let cell parameters free and monitor the VOLUME,
then use the formula above.
Or - probably better - you extract elastic constants,
checking back their definition in hcp, as indepent parameters.
Because they may behave differently; the bulk modulus is just
some averaged combination of them, and won't tell you much,
for a serious test.

>> is it enough to use md.variable cell? how will i know that cell geometry
>> remains hexagonal?
>
>  Of course you shouldn't use variable cell: for each \delta V above the
>  cell is kept fixed, otherwise the pressure will relax to zero.

Sorry again,
I thing, it is ESSENTIAL to use variable cell, for the reason given above.

>> should i fix the atomic positions?
>
>  No. The hcp 2-atoms cell possesses inner degrees of freedom,

How's that??? I always thought there is something like
(0 0 0) and (1/3  2/3  1/2), no internal coordinates.
These relative coordinates should not change
(they will in fact, slightly, due to lack of symmetry constraint
in Siesta).

> most likely the symmetry will be broken.
>
>> Is it all right if i test the pseudo for fcc structure instead and
>> compare
>> the results with available theoretical results.

And what you want to compare with otherwise, the experiment?
Is there some, for technetium? Just wondering.

>  Who knows. Any test is just a test, it should be helpful anyway.

Hmm. I think, having a good all-electron code at hand to compare,
you can, to begin with, construct any distortion you want,
do the same with Siesta, compare one curve against the other
and already get an idea.
For this you don't even need to know what the bulk modulus is.

Best regards

Andrei