Re: [SIESTA-L] Regarding SCF convergence issue during optimization

[Andrei Postnikov]

Hi, 
you provide too few information. Of course just increasing the number of 
iterations won't bring you further 
if there is no convergence. General advises might be as follows: 

- Check that you are using diagonalisation and NOT order-N scheme before you 
fix the things and 
"know what you are doing". With 235 atoms the chances are, you are in order-N 
regime by default, 
and things go astray. 

- Check your system visually (using visualisation tools from your working XV 
file) 
whether its shape is as expected; a huge number of problems result from an 
error in defining the structure. 

- Check whether your electronic structure makes sense (the main bands are more 
or less where they are expected) 
at very first iterations, before the calculations start to markedly diverge. Do 
not just look at the convergence, but at other things 
in the output file. Make sure that the attribution of valence states in the 
pseudopotential and in the basis matches. 
Some inconsistencies may cry aloud in the output file, but you need to pay 
attention... 

- Reduce the mixing parameter. Reduce it drastically if needed (0.001 is not 
uncommon in some cases). 
Start converging with large broadening (electronic temperature), say 600 K, and 
then reduce it gradually as you enter 
into a convergence regime. 

The actual behaviour may depend on whether your system is isolated or periodic, 
insulator or metal. 
If your case is periodic and metallic, you may need many k-points as well (in 
the directions of periodicity, of course). 

Good luck 

Andrei Postnikov 

- Le 19 Avr 24, à 4:16, Jyotirmoy Deb  a écrit : 

> Dear Sir/Madam,
> I want to optimize a heterostructure consisting of 235 atoms. But I am facing 
> an
> error every time "SCF did not converge in the maximum number of steps". I have
> gradually increased the "MaxSCFIterations" however, I have not been able to
> solve this issue. Please help me in this regard.

> Thanking you
> with regards
> Jyotirmoy

> Dr. Jyotirmoy Deb
> C/O Dr. G. Narahari Sastry & Dr. Lakshi Saikia
> DST-SERB National Post-Doctoral Fellow
> Advanced Computation and Data Sciences Division (ACDSD)
> CSIR-North East Institute of Science & Technology, Jorhat-785006, Assam, 
> India.
> Ph. No: +917002140643/ +919435589869
> Email: [ mailto:deb.jyotirmo...@gmail.com | deb.jyotirmo...@gmail.com ]
> Webpage: [
> https://urldefense.com/v3/__https://sites.google.com/view/jyotirmoydftphy__;!!D9dNQwwGXtA!QXMXS7L8cCba3hBh8KgcbwyC9C-E5X98KJTAXUguv0CMpKRil5P5cMWK5ziJhwU4-c14ht8eejcbtj_czqW_iUBo$
> | 
> https://urldefense.com/v3/__https://sites.google.com/view/jyotirmoydftphy__;!!D9dNQwwGXtA!TlNGLjwuIHLvaLUEXMAzZneDEOzZtFJFMRfYAnBQADgcqVNh1tJuG4DMBR8sWhbOVnrRtR-96vIssOhkbvkqvQKGf5LpgmFIQA$
>   ]

> --
> SIESTA is supported by the Spanish Research Agency (AEI) and by the European
> H2020 MaX Centre of Excellence 
> (https://urldefense.com/v3/__http://www.max-centre.eu/__;!!D9dNQwwGXtA!TlNGLjwuIHLvaLUEXMAzZneDEOzZtFJFMRfYAnBQADgcqVNh1tJuG4DMBR8sWhbOVnrRtR-96vIssOhkbvkqvQKGf5IASVdD0A$
>  )

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Problem using Vibra utility

[Andrei Postnikov]

Dear Diego,

as your case is q=0 only, you can try your luck
with my shortcut version of vibra, to be compiled from the attached zip.
It does not need an .fdf file, just .XV and .FC 
(.XV serves just to identify the atoms; exact coordinates are not important).
Tell me if you'd encounter any difficulties.

Best regards

Andrei Postnikov 


- Le 12 Sep 23, à 12:16, Diego Lopez Alcala diego.lo...@uv.es a écrit :

> Dear Siesta users,
> 
> I have been trying to compute the modes of vibrations of a molecule adsorbed 
> on
> a semiconducting monolayer, but I am having some problems. First I run this
> input:
> 
> # General System descriptors
> 
> SystemName vibra-2   # Descriptive name of the system
> SystemLabel   vibra-2# Short name for naming files
> 
> NumberOfAtoms   164   # Number of atoms
> NumberOfSpecies 5   # Number of species
> 
> md.typeofrun fc
> 
> MD.FCFirst 151
> MD.FCLast 164
> 
> AtomicCoordinatesFormat NotScaledCartesianAng
> 
> MM.UnitsDistance Ang  # what this program prints out DO NOT CHANGE
> MM.UnitsEnergyeV  # what this program prints out DO NOT CHANGE
> MM.Grimme.S6 0.75 # Grimme-paper for PBE (correct for your functional)
> MM.Grimme.D 20.   # Grimme-paper (correct for your functional)
> %block MM.Potentials
>  1   1 Grimme111.94  3.124 # Cr, 10.1002/jcc.20495
>  1   2 Grimme 80.39  3.245 # Cr / S
>  1   3 Grimme120.28  3.311 # Cr / Br
>  1   4 Grimme 45.06  3.014 # Cr / C
>  1   5 Grimme 12.74  2.563 # Cr / H
>  2   2 Grimme 57.73  3.366 # S, 10.1002/jcc.20495
>  2   3 Grimme 86.38  3.432 # S / Br
>  2   4 Grimme 32.36  3.135 # S / C
>  2   5 Grimme  9.15  2.684 # S / H
>  3   3 Grimme129.24  3.498 # Br, 10.1002/jcc.20495
>  3   4 Grimme 48.42  3.201 # Br / C
>  3   5 Grimme 13.69  2.750 # Br / H
>  4   4 Grimme 18.14  2.904 # C, 10.1002/jcc.20495
>  4   5 Grimme  5.13  2.453 # C / H
>  5   5 Grimme  1.45  2.002 # H, 10.1002/jcc.20495
> %endblock MM.Potentials
> 
> %block Chemical_Species_Label
>  1   24 Cr
>  2   16 S
>  3   35 Br
>  46 C
>  51 H
> %endblock Chemical_Species_Label
> 
> PAO.BasisSize  SZ
> 
> DFTU.ProjectorGenerationMethod 1
> DFTU.PotentialShift .true.
> 
> %block DFTU.Proj # Define DFTU projectors
> Cr 1 # Label, l_shells
> n=3 2  # n (opt if not using semicore levels),l,Softconf(opt)
> 4.00 0.0 # U(eV), J(eV) for this shell
> 0.0 # rc (Bohr), \omega(Bohr) (Fermi cutoff function)
> %endblock DFTU.Proj
> 
> # Lattice, coordinates, k-sampling
> 
> AtomicCoorFormatOut Ang
> 
> %block AtomicCoordinatesAndAtomicSpecies
> 4.44717269  3.68308271  0.75902034  3   79.904
> .
> . (Rest of atomic coordinates)
> .
> 5.10327585  14.06972694 8.19442056  5   1.007
> %endblock AtomicCoordinatesAndAtomicSpecies
> 
> LatticeConstant 1.0 Ang
> 
> %block LatticeVectors
>   17.8287990.000.00
>0.00   24.3147200.00
>0.000.00   25.236237
> %endblock LatticeVectors
> 
> %block kgrid_Monkhorst_Pack
> 1   0   0  0.0
> 0   1   0  0.0
> 0   0   1  0.0
> %endblock kgrid_Monkhorst_Pack
> 
> # DFT, Grid, SCF
> 
> XC.functional   GGA # Exchange-correlation functional type
> XC.authors  PBE # Particular parametrization of xc func
> SpinPolarized   .true.  # Spin unpolarized calculation
> MeshCutoff  400. Ry # Equivalent planewave cutoff for the grid
> MaxSCFIterations300 # Maximum number of SCF iterations per 
> step
> SCF.DM.Converge true
> SCF.H.Converge  true
> 
> # Eigenvalue problem: order-N or diagonalization
> 
> SolutionMethod  diagon  # OrderN or Diagon
> ElectronicTemperature   300 K# Temp. for Fermi smearing
> 
> # Output options
> 
> WriteCoorInitial
> WriteCoorStep
> WriteForces
> WriteMullikenPop1 # Write Mulliken Population Analysis
> WriteCoorXmol   .false.
> WriteMDCoorXmol .false.
> WriteMDhistory  .false.
> WriteCoorXmol   .false.
> 
> Once it is done I add the following lines to the input (as suggested  here:
> https://urldefense.com/v3/__https://docs.siesta-project.org/projects/siesta/en/latest/tutorials/basic/vibrational-properties/Benzene/index.html__;!!D9dNQwwGXtA!UewEzyC8ZoPpI3i_PmoF3Vbbcs86V_CEf3F-hofbAGex6SS1dlsBvcgSPg7HdC2VlFmiZZGfJEXGqNlpcyY$
> ) and I run the vibra utility, obtaining this error message. The
> SystemLabel.bands is

Re: [SIESTA-L] Radial part of basis orbital

[Andrei Postnikov]

Isn't this information explicitly included in .ion files - 
after a general header, one basis function after the other? 
You just need to cut whichever you want and feed it to a plotting program... 
Best regards, 
A.P. 

... 

# PAOs:__ 

0 2 1 0 2.00 #orbital l, n, z, is_polarized, population 

500 0.132468670984E-01 6.61018668210 # npts, delta, cutoff 

0.000 0.287172515330 

0.132468670984E-01 0.287187021634 

0.264937341968E-01 0.287232752566 

0.397406012952E-01 0.287308980728 

0.529874683936E-01 0.287415720365 

... 

- Le 15 Fév 23, à 15:34, Emilio Artacho  a écrit : 

> Check the WriteIonPlotFiles option

> E

>> On 14 Feb 2023, at 06:17, Francisco Garcia < [ 
>> mailto:garcia.ff@gmail.com |
>> garcia.ff@gmail.com ] > wrote:
>> Dear Users,

>> I would like to know how to obtain the data for the radial part of a given 
>> basis
>> orbital for plotting. E.g. the radial part of Al 3s and Al 3p for SZ basis.

>> Thank you.

>> --
>> SIESTA is supported by the Spanish Research Agency (AEI) and by the European
>> H2020 MaX Centre of Excellence ( [
>> https://urldefense.com/v3/__http://www.max-centre.eu/__;!!D9dNQwwGXtA!S0apBc_vxFM3XsCbZ2x_Ek6V_jNWS_G7uxbf6fayBqjsXL1bU6GTY35UzarXF3pSaMhFMFRLKKd94oJ0dFV9$
>> | 
>> https://urldefense.com/v3/__http://www.max-centre.eu/__;!!D9dNQwwGXtA!W61yKVCs4FflUUO_VbPHAuO0S3HZ8Xzc17Dr0DE6UJW862iIwupaYKBiGgBCBcoIWd0sdJ20S4blr36ZkIE3e-aI_Z3qKyHZtQ$
>>   ] )

> --
> Emilio Artacho

> Theory Group, Nanogune, 20018 San Sebastian, Spain, and
> Theory of Condensed Matter, Department of Physics,
> Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, UK

> --
> SIESTA is supported by the Spanish Research Agency (AEI) and by the European
> H2020 MaX Centre of Excellence 
> (https://urldefense.com/v3/__http://www.max-centre.eu/__;!!D9dNQwwGXtA!W61yKVCs4FflUUO_VbPHAuO0S3HZ8Xzc17Dr0DE6UJW862iIwupaYKBiGgBCBcoIWd0sdJ20S4blr36ZkIE3e-aI_Z3qKyHZtQ$
>  )

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Restarting calculations with different spin orientations

[Andrei Postnikov]

Dear Andres Tellez, 
my DMtune tool was written having similar tasks in mind; 
please feel free to download it from 
[ 
https://urldefense.com/v3/__https://www.home.uni-osnabrueck.de/apostnik/download.html__;!!D9dNQwwGXtA!UOZi7W_dTZOE_UPIXqgD4MdGacy1Fx8Kvs3HymrL_qKWalh8UtdNSwuKc4JOVzfzoPLVLIAZ5uHSbQoT-gAurVYvIzKoJQhPng$
  | 
https://urldefense.com/v3/__https://www.home.uni-osnabrueck.de/apostnik/download.html__;!!D9dNQwwGXtA!UOZi7W_dTZOE_UPIXqgD4MdGacy1Fx8Kvs3HymrL_qKWalh8UtdNSwuKc4JOVzfzoPLVLIAZ5uHSbQoT-gAurVYvIzKoJQhPng$
  ] 
and test whether it still works and applies to your case. 
There is no special documentation, but minimal explanations are included 
at the beginning of the source file dmtune.f. 
In case of problems it should not be difficult to modify, or pass me a word 
(with an example), I'll try to update. 

Best regards 

Andrei Postnikov 

- Le 31 Jan 23, à 18:14, Andres Tellez Mora  a écrit : 

> Dear Siesta users and developers.

> I am running calculations of the same structure with different spin
> orientations. Since I am performing spin-orbit calculations with DFT+U, they
> can take a decent amount of time to converge. Hence, I was trying to use the
> .DM file of one of the calculations to start the others; however, the
> DM.InitSpin block information is overridden and the calculation starts using
> the same spin orientation from the converged .DM file. This even happens when
> using a .DM file from a non-polarized calculation, which gives a starting
> magnetization of 0.0. Is it possible to use the information of a .DM file and
> start with a given spin orientation simultaneously? If this is not possible,
> what else could I do to improve the convergence? I appreciate any help you can
> provide.

> Best Regards,
> Andres Tellez.

> --
> SIESTA is supported by the Spanish Research Agency (AEI) and by the European
> H2020 MaX Centre of Excellence 
> (https://urldefense.com/v3/__http://www.max-centre.eu/__;!!D9dNQwwGXtA!UOZi7W_dTZOE_UPIXqgD4MdGacy1Fx8Kvs3HymrL_qKWalh8UtdNSwuKc4JOVzfzoPLVLIAZ5uHSbQoT-gAurVYvIzK9M_8LNA$
>  )

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Need help with Phonon Dispersion Band Lines in Vibra

[Andrei Postnikov]

Dear Francisco, 
so your AFM structure is of CuAu type. This is fine but of course 
this is not the only AFM structure possible (and I don't know whether it is 
realistic at all, but this of course depends on your objectives). 
Now, if you want the q-path to be the same in your two settings, 
you should consider that the second one is rotated by 45 degrees. 
That is, if you choose the Gamma->X direction || [010] in the first setting 
it must be || [110] in the second setting, with the lattice vectors you use. 
Otherwise define the lattice vectors as [1/2 -1/2 0], [1/2 1/2 0], [0 0 1] 
with the same lattice parameter as in the first setting 
and enjoy the same coordinates of q points (cartesian, in terms of pi/a) 
in both settings. 

Best regards 

Andrei 

- Le 31 Déc 22, à 0:24, garcia ff 000  a écrit : 

> Dear Prof. Postnikov,

> Many thanks and appreciation for your response. I believe I found a solution 
> to
> my problem but I want to run it by you.

> First, an FCC cell with 2 unique atoms is equivalent to a tetragonal cell 
> (this
> is the smallest unit cell to model antiferromagnetism).

> Using the website [ 
> https://urldefense.com/v3/__https://www.materialscloud.org/work/tools/seekpath__;!!D9dNQwwGXtA!VCzh9W4S1t5nhfeK_65w_ZsZpJauei8vdCoYcoysbxbXQ6kbxNBuSTzR-LciHx145nkwK_JGfqplTyD_aZeQ8icIeWJtZ4jCPg$
>   |
> https://urldefense.com/v3/__https://www.materialscloud.org/work/tools/seekpath__;!!D9dNQwwGXtA!VCzh9W4S1t5nhfeK_65w_ZsZpJauei8vdCoYcoysbxbXQ6kbxNBuSTzR-LciHx145nkwK_JGfqplTyD_aZeQ8icIeWJtZ4jCPg$
>   ] , the high symmetry points
> in the Brillouin zone are as follows (each set of points is scaled by the
> corresponding pi/a):

> Standard FCC primitive cell: Gamma (0,0,0), X(0,2,0), K(1.5,1.5,0), W(1,2,0),
> L(1,1,1)

> 2-atom tetragonal cell: Gamma(0,0,0), X(0,1,0), M(1,1,0), R(0,1,0.707107),
> A(1,1,0.707107), Z(0,0,0.707107).

> With this information, I believe the two Vibra inputs below, one for the
> primitive FCC cell and the other for 2-atom tetragonal cell, are formally
> equivalent (the last two k-points in each case, i.e. L and M, is what I'm a 
> bit
> unsure about).

> Thank you very much for your kindness & happy holidays.

> (A) Primitive FCC cell:

> NumberOfAtoms 1

> #Lattice parameters
> LatticeConstant 3.47 Ang
> %block LatticeVectors
> 0.50 0.50 0.00
> 0.50 0.00 0.50
> 0.00 0.50 0.50
> %endblock LatticeVectors

> #Atomic positions
> AtomicCoordinatesFormat Fractional
> %block AtomicCoordinatesAndAtomicSpecies
> 0.00 0.00 0.00 1 54.938
> %endblock AtomicCoordinatesAndAtomicSpecies

> #High symmetry Brillouin zones points scaled by pi/a: Gamma (0,0,0), X(0,2,0),
> K(1.5,1.5,0), W(1,2,0), L(1,1,1)

> BandLinesScale pi/a
> %block BandLines
> 1 0.000 0.000 0.000 \Gamma
> 30 0.000 2.000 0.000 X
> 30 2.000 2.000 2.000 \Gamma
> 30 1.000 1.000 1.000 L
> %endblock BandLines

> (B) 2-atom tetragonal cell to model antiferromagnetism (this is double the
> volume of the FCC primitive cell)

> NumberOfAtoms 2

> #Lattice parameters
> LatticeConstant 2.453660531 Ang #[this is the FCC lattice constant divided by
> sqrt(2)]
> %block LatticeVectors
> 1.00 0.00 0.00
> 0.00 1.00 0.00
> 0.00 0.00 1.414213562
> %endblock LatticeVectors

> #Atomic positions
> AtomicCoordinatesFormat Fractional
> %block AtomicCoordinatesAndAtomicSpecies
> 0.00 0.00 0.00 1 54.938
> 0.50 0.50 0.50 1 54.938
> %endblock AtomicCoordinatesAndAtomicSpecies

> #High symmetry Brillouin zones points scaled by pi/a: Gamma(0,0,0), X(0,1,0),
> M(1,1,0), R(0,1,0.707107), A(1,1,0.707107), Z(0,0,0.707107)

> BandLinesScale pi/a
> %block BandLines
> 1 0.000 0.000 0.000 \Gamma
> 30 0.000 1.000 0.000 X
> 30 1.000 1.000 1.000 \Gamma
> 30 2.000 2.000 2.000 M
> %endblock BandLines

> On Thu, Dec 29, 2022 at 3:34 PM Andrei Postnikov < [
> mailto:andrei.postni...@univ-lorraine.fr | andrei.postni...@univ-lorraine.fr ]
> > wrote:

>> Dear Francisco,
>> it is difficult to give a useful advice on the basis of very limited 
>> information
>> you provide,
>> but my impression is that your problems are not obviously related with Vibra.
>> Some questions:
>> 1. What (magnetic) structure are you modelling? How comes you have four atoms
>> per AFM unit cell?
>> Can there be two?
>> 2. Is electronic structure (and band dispersions) correct, prior to any 
>> phonons?
>> 3. What means "incorrect phonon dispersion"? Do you have problems with
>> crystallography /
>> choosing the q-path, or is your calculation basically wrong?
>> 4. With 4 atoms as you use it so far, the Gamma phonon calculation would 
>>

Re: [SIESTA-L] Need help with Phonon Dispersion Band Lines in Vibra

[Andrei Postnikov]

Dear Francisco, 
it is difficult to give a useful advice on the basis of very limited 
information you provide, 
but my impression is that your problems are not obviously related with Vibra. 
Some questions: 
1. What (magnetic) structure are you modelling? How comes you have four atoms 
per AFM unit cell? 
Can there be two? 
2. Is electronic structure (and band dispersions) correct, prior to any 
phonons? 
3. What means "incorrect phonon dispersion"? Do you have problems with 
crystallography / 
choosing the q-path, or is your calculation basically wrong? 
4. With 4 atoms as you use it so far, the Gamma phonon calculation would yield 
9 modes, which would map genuine zone-center and zone-boundary modes. 
Do they come out reasonably? 

To your problem: 
"B asically I want to alter the band lines in input 2 so that they are 
equivalent to the band lines in input 1" - 
you have 
BandLinesScale pi/a 
in both inputs, the same lattice parameter, and the same definition of path. 
So if everything is correctly read, you must get the same Cartesian q-path in 
both cases. 
Either this is not so and there is something wrong with the input, 
or the paths are identical but your problem is elsewhere. 

Best regards 

Andrei 

- Le 29 Déc 22, à 0:40, garcia ff 000  a écrit : 

> Dear Users,

> I have appended 2 Vibra inputs below for computing the phonon dispersion for 
> FCC
> Mn.

> Input 1 works fine as it gives the expected band shapes for the dispersion 
> (but
> the frequencies are off). The main issue with input 1 is that it is not
> suitable for antiferromagnetic calculations since there is only one Mn atom in
> the primitive cell.

> This led me to consider input 2, which has 4 atoms in the unit cell and can be
> used for antiferromagnetic calculations. The issue with input 2 is that the
> bandlines yield an incorrect phonon dispersion. This is what I need your help
> on. Basically I want to alter the band lines in input 2 so that they are
> equivalent to the band lines in input 1.

> Any assistance with this, especially from the Vibra authors, would be greatly
> appreciated.

> Thank you very much for your kind assistance and God Bless!

> Francisco

> #INPUT 1 (1 atom in the FCC primitive cell; 125 atoms in Supercell)
> SystemName fccMn_1
> SystemLabel fccMn_1
> NumberOfAtoms 1
> LatticeConstant 3.47 Ang
> %block LatticeVectors
> 0.50 0.50 0.00
> 0.50 0.00 0.50
> 0.00 0.50 0.50
> %endblock LatticeVectors

> AtomicCoordinatesFormat Fractional
> %block AtomicCoordinatesAndAtomicSpecies
> 0.00 0.00 0.00 1 54.938
> %endblock AtomicCoordinatesAndAtomicSpecies

> SuperCell_1 2
> SuperCell_2 2
> SuperCell_3 2

> AtomicDispl 0.04 Bohr

> BandLinesScale pi/a
> %block BandLines
> 1 0.000 0.000 0.000 \Gamma
> 30 2.000 0.000 0.000 X
> 30 2.000 2.000 2.000 \Gamma
> 30 1.000 1.000 1.000 L
> %endblock BandLines

> Eigenvectors True

> #INPUT 2 (4 atoms in the FCC conventional cell; 108 atoms in Supercell)
> SystemName fccMn_4
> SystemLabel fccMn_4
> NumberOfAtoms 4
> LatticeConstant 3.47 Ang
> %block LatticeVectors
> 1.00 0.00 0.00
> 0.00 1.00 0.00
> 0.00 0.00 1.00
> %endblock LatticeVectors

> AtomicCoordinatesFormat Fractional
> %block AtomicCoordinatesAndAtomicSpecies
> 0.00 0.00 0.00 1 54.938
> 0.50 0.50 0.00 1 54.938
> 0.50 0.00 0.50 1 54.938
> 0.00 0.50 0.50 1 54.938
> %endblock AtomicCoordinatesAndAtomicSpecies

> SuperCell_1 1
> SuperCell_2 1
> SuperCell_3 1

> AtomicDispl 0.04 Bohr

> BandLinesScale pi/a
> # The band lines below are incorrect.
> %block BandLines
> 1 0.000 0.000 0.000 \Gamma
> 30 2.000 0.000 0.000 X
> 30 2.000 2.000 2.000 \Gamma
> 30 1.000 1.000 1.000 L
> %endblock BandLines

> Eigenvectors True

> --
> SIESTA is supported by the Spanish Research Agency (AEI) and by the European
> H2020 MaX Centre of Excellence 
> (https://urldefense.com/v3/__http://www.max-centre.eu/__;!!D9dNQwwGXtA!T-wl-ZvX-LX5xZC7QdfhJRIJ8Pmxo5HofWGt13XzKiGWhx9VgP3MXmjkoKcw2oYTy4STEEQIyWW5lU0aV4mzNGNhq7rtk7d9KA$
>  )

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re:[SIESTA-L]

[Andrei Postnikov]

Dear Amal Yassin: 
Your unit cell vectors are all zeros; 
your atom belongs to specie Nr zero (must be 1), 
and there seems to be a problem with its coordinate. 
All this is evident from your output file if you read it beyond the first line. 
The absence of XV file is just a warning; it's no problem 
if the structure is correctly defined in the input file. 

Good luck, 

Andrei Postnikov 

- Le 18 Aoû 22, à 22:18, Amal Yassin  a écrit : 

> Hello
> When i run the simulation of Carbone , i existe this error .
> Any Idea about this error , and how can resolve it ?!
> Please Help..

> Best regard

> siesta: WARNING: XV file not found

> siesta: Atomic coordinates (Bohr) and species
> siesta: 0.0** 3.77945 0 1

> siesta: Automatic unit cell vectors (Ang):
> siesta: 0.00 0.00 0.00
> siesta: 0.00 0.00 0.00
> siesta: 0.00 0.00 0.00
> rcut: Wrong species 0. Have 1
> Stopping Program from Node: 0
> #0 0x7f0bf0cdbd21 in ???
> #1 0x7f0bf0cdc7ad in ???
> #2 0x7f0bf0f3038c in ???
> #3 0x55c259359aca in ???
> #4 0x55c259098e83 in ???
> #5 0x55c25901ed25 in ???
> #6 0x55c25904bd0b in ???
> #7 0x55c25904a277 in ???
> #8 0x55c258f8f74c in ???
> #9 0x7f0bf0957082 in __libc_start_main
> at ../csu/libc-start.c:308
> #10 0x55c258f8f81d in ???
> #11 0x in ???

> --
> SIESTA is supported by the Spanish Research Agency (AEI) and by the European
> H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Forces on crystalline atoms are not zero

[Andrei Postnikov]

Dear Francisco:as a general statement, your expectation is justified. However, 
unexpected things happen.Having posted an (input) / output file would help to 
pinpoint a problem. Without any additional information, I'd guess that your 
lattice vectors are not exactly fcc, or the atoms are not exactly at symmetric 
positions. (E.g., you think that you scalethe atom coordinates with lattice 
parameter, but in fact you don't).Best regardsAndrei Postnikov
- Francisco Garcia  a écrit :
>Dear Users
I performed a series of single point energy calculations on an FCC crystal by 
varying the lattice constant. I'm trying to generate data points to optimize 
the lattice constant

I was expecting the force on each atom to zero by virtue of crystalline 
symmetry (no internal degree of freedom). But to my surprise, the forces are 
not zero. Increasing the Mesh cutoff from 300 Ry to higher values didn't help 
either. How is it possible for the atomic forces in a crystal to be non-zero?



-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] phdos

[Andrei Postnikov]

Right.Sorry if this was not clear enough from the documentation:my "phdos" is 
intended to be used on Gamma phonon calculated on a (probably large) cell,that 
is, for a case when the number of atoms is the .XV file is the same asthe 
number of atoms for which eigenvectors are calculated.The idea is, to broaden a 
discrete spectrum of (presumably many) linesof Gamma-vibration frequencies, not 
to integrate the density of modes fromthe phonon dispersions over the Brillouin 
zone.In your case, the .XV file contains 125 atoms (5x5c1 supercell à 5 
atoms)but the eigenvectors are calculated for 5 atoms only (and this over a 
path of many q-points).Ghe "phdos" identifies 125 atoms in the .XV file and 
attempts to read the eigenvectors for all of thembut finds the data for 5 atoms 
only, hence the error message and termination.So such behaviour is as 
expected.In order to solve exactly your problem I started to work some time ago 
on an auxiliary code which expandsa "partial" FC file into a "full" FC file 
(which would result from displacing all atoms in the supercelland not only 
"prototype" ones in the original single cell). It was not extensively testedand 
does not immediately work on our example although in principle it should; I'll 
check what is the problem and let you know.Best regardsAndrei Postnikov
- AK- HYDRA  a écrit :
>Hello sir/ma'amHere, I got same issue every time i.e., Error reading real part 
>of eigenvector  1  iat=6while using 'phdos' available on github within 
>'Lattice_dynamics'.I am going to attach my .vector, .FC, *.fdf, and .XV file, 
>need suggestion please 
thank you 

-- 


-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] strange force in phonon calculation

[Andrei Postnikov]

Dear Ziba,it seems that you are randomly trying different options.Try to be 
more systematic.MeshCutoff of 870 Ry seems to be an overkill -but certainly not 
an error; if your tests show that you need as much as that, so be it.Did you 
have a look at your structures with visualisation tools - is everything as it 
should be, no obvious structure error?If your imaginary frequency is genuine 
and reproducible,I'd look at the eigenvector of this mode, give a small 
displacementto atoms along this eigenvector, and make the relaxation again. The 
total energy ought to be lower than in your initial starting point;this is the 
idea of an imaginary frequency: you displace the atomswithin the corresponding 
mode, and the energy goes down.Hence you were not at the equilibrium; try 
better.No way around.Best regardsAndrei Postnikov
- Ziba Torkashvand  a écrit :
>
Thank you Professor Postnikov,Here I will provide two examples of my workI have 
a hexagonal structure containing 18 atoms1. I have substituted one atom with a 
gaseous one and relaxed it up to 0.000794 eV/A and 0.00458798  kBarthen I 
calculated the phonon and fortunately, there is no negative frequency while the 
force for the first step is approximated as 0.008174 eV/A
2. I have substituted two atoms with two gaseous ones and relaxed it up to 
0.000549 eV/A and 0.00290065  kBarthen I calculated the phonon and 
unfortunately, there are frequencies as -206.30292747012965 while the force for 
the first step is approximated as 0.015862 eV/A
For both cases, I have used same inputs asFor relaxation:MeshCutoff870 Ry
DM.MixingWeight 0.02
DM.NumberPulay   3
%block kgrid_Monkhorst_Pack
6   0   0   0.0
0   6   0   0.0
0   0   1   0.0
%endblock kgrid_Monkhorst_Pack

For phonon: MeshCutoff750 Ry
DM.MixingWeight 0.02
DM.NumberPulay   8
%block kgrid_Monkhorst_Pack
3   0   0   0.0
0   3   0   0.0
0   0   1   0.0
%endblock kgrid_Monkhorst_Pack

I repeated phonon calculation for the second case with %block 
kgrid_Monkhorst_Pack
6   0   0   0.0
0   6   0   0.0
0   0   1   0.0
%endblock kgrid_Monkhorst_Pack
and the resulting force for the first step is 0.013928 eV/A
I'm confused about what is happening in these two similar examples.

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] strange force in phonon calculation

[Andrei Postnikov]

Dear Ziba,probably you'll need to explain better what's your problem.> Once I'm 
relaxing my structure up to 0.0004 eV/Ang after the first MD step of phonon 
calculation > without any displacement the maximum force will be 0.002 eV/Ang 
"The first step without any displacement" is just calculation of your relaxed 
structure.If it is as you say then either it was not fully relaxed, or you 
changed some parametersbetween the two calculations. (A guess: either the 
coordinates are not those from the relaxed case,or the density matrix taken was 
not converged).> but still there is no negative frequency.- what you mean, 
"still"? And how do you know "there is no negative frequency" if you've just 
done "the first step without any displacement"?> This time I've relaxed my 
structure up to 0.0002 eV/Ang but in the first step I'm receiving 0.015 eV/Ang- 
see above> which is leading to the negative frequencies.- for what I know, 
negative frequencies are are hint that the starting point for a phonon 
calculationwas not at the local minimum, i.e., the relaxation was not good. How 
negative are your frequencies?> In the successful phonon calculations I just 
lowered the number of kpoint in the Monkhorst pack for phonon calculation > but 
this time even the kpoints are the same as relaxation input.- reducing the 
k-mesh is certainly not a working recipe for getting phonons right; you've just 
had luck.Best regardsAndrei Postnikov
- Ziba Torkashvand  a écrit :
>
Hello everyone,I'm trying to do phonon calculation for some hexagonal 
structures.Once I'm relaxing my structure up to 0.0004 eV/Ang after the first 
MD step of phonon calculation without any displacement the maximum force will 
be 0.002 eV/Ang but still there is no negative frequency.This time I've relaxed 
my structure up to 0.0002 eV/Ang but in the first step I'm receiving 0.015 
eV/Ang which is leading to the negative frequencies.In the successful phonon 
calculations I just lowered the number of kpoint in the Monkhorst pack for 
phonon calculation but this time even the kpoints are the same as relaxation 
input.
Any help would be appreciated
ThanksZiba

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Query regarding charge calculation using SCF charge density difference

[Andrei Postnikov]

Dear Harkishan,I think Macroave (in the Util/ of the Siesta distribution)might 
be able to do something similar to what you need.Otherwise, I've a simple 
script (attached) intended to integrate(in fact, simply to sum up) the charge 
density given on the grid over spherescentered on atoms. The idea was to get a 
fair comparison with chargesestimated in other ("muffin-tin") codes. This 
worked fairy well even if not very fast.You can modify it to integrate over a 
user-defined volume.The algorithm simply runs over all points of the grid to 
checkwhether they fall within the defined shape.Hopefully this may help,best 
regardsAndrei Postnikov
- Harkishan Dua  a écrit :
>I have a bilayer system which is supposed to act as a nanoscale capacitor. I 
>have performed the calculations using siesta where I have calculated the 
>charge stored in the capacitor from the dipole moment. But after communicating 
>the paper, the reviewers suggested getting the charge stored in the capacitor 
>using SCF charge density difference. Please guide me on how to get the charge 
>stored in the capacitor using SCF charge density difference in siesta. I would 
>be really grateful for any help in this regard.
Eagerly waiting for your reply
 
Regards
Harkishan Dua
-- 
Harkishan DuaPhD StudentDepartment of PhysicsAssam University, Silchar.C
Cgrdint,  a script to integrate 3-dim grid function
C(i.e. LDOS, RHO, DRHO, etc.) written in SIESTA by subr. iorho
C over a region of choice (a sphere centered at a given atom).
C
C   !!!   IMPORTANT  -- !!!
C   compile this code with the same compiler switches as Siesta,
C in what regards using single/double precision,
C   otherwise reading the data from unformatted files 
C  written by iorho.F  might be spoiled.  
CDon't say you haven't been warned.
C!!!
C
C Written by Andrei Postnikov, Mar 2007  
C postni...@univ-metz.fr
C
  program grdint
  implicit none
  integer ii1,ii2,io1
  parameter (ii1=11,ii2=12,io1=14)
  integer mesh(3),nspin,is,ii,jj,iat,nat,nz,iix,iiy,iiz,
 .ind,mn,ialloc,ishift,jshift,kshift,limit,nrat
  parameter (limit=1)  !  tried translations along each lattice vector
  character inpfil*60,outfil*60,syslab*30,suffix*6,owrite*1
  logical filexist
  double precision b2ang,cc_bohr(3,3),diff,coort(3),cell(3,3),
 . dummy,modu,rmesh(3),small,volume,srad,srad2,
 . dist2,charge,spin
  parameter (b2ang=0.529177)   !  Bohr to Angstroem
  parameter (small=1.0d-5)  
  integer, allocatable :: ityp(:),iz(:)
  double precision, allocatable :: coor_bohr(:,:)
  real, allocatable :: func(:,:,:,:) !  NB! single precision, as in iorho.F
  real  sum(8)   !  NB! single precision, as in iorho.F
C
C string manipulation functions in Fortran used below:
C len_trim(string): returns the length of string 
C   without trailing blank characters,
C trim(string): returns the string with railing blanks removed
C char(integer)   : returns the character in the specified position
C   of computer's ASCII table, i.e. char(49)=1
  
  write (6,701,advance="no")
  701 format(" Specify  SystemLabel (or 'siesta' if none): ")
  read (5,*) syslab
C --  handle XV file: 
  inpfil = trim(syslab)//'.XV'
  open (ii1,file=inpfil,form='formatted',status='old',err=801)
C --  check number of atoms in the XV file, for allocation:
  do ii=1,3
read (ii1,'(3x,3f18.9)',end=801,err=801) (dummy,jj=1,3)
  enddo
  read (ii1,*,end=802,err=802) nat
  allocate (ityp(nat))
  allocate (iz(nat))
  allocate (coor_bohr(1:3,1:nat),STAT=ialloc)
  if (ialloc.ne.0) then
 write (6,*) ' Fails to allocate space for ',nat,' atoms.'
 stop
  endif
  rewind (ii1)
C --  read in translation vectors: 
  do ii=1,3
read  (ii1,'(3x,3f18.9)')  (cc_bohr(jj,ii),jj=1,3)
  enddo
  read (ii1,*) nat
  do iat=1,nat
read (ii1,'(i3,i6,3f18.9)',end=803,err=803) 
 .ityp(iat),iz(iat),coor_bohr(:,iat)
  enddo
  close (ii1)
C
C --- Look for grid data files to include:
  103 write (6,706,advance="no")
  read (5,*) suffix
  inpfil = trim(syslab)//'.'//trim(suffix)
  open (ii2,file=inpfil,form='unformatted',status='old',err=806)
  write (6,*) 'Found and opened: ',trim(inpfil)
  read (ii2,err=807) cell
C check if cell vectors match those from .XV, just for the case...
  diff = 0.0
  do ii=1,3
  do jj=1,3
diff = diff + (cell(ii,jj)-cc_bohr(ii,jj))**2
  enddo
  enddo
  if (diff.gt.small) then
write (6,*) ' cell vectors from XV: '
write (6,'(3f12.6)') (cc_bohr(ii,

Re: [SIESTA-L] Band structure for BiFeO3

[Andrei Postnikov]

Dear ionutghitiu,your message is a bit confusing. You say your band structures 
are like in other calculations,the magnetic moments are like in the literature 
- what is the problem, then?The list of different parameters you tried is 
impressive but doesn't give a clue.You have large Fe magnetic moments but you 
see no spin splitting in (total?) DOS and bands -are you looking at an 
antiferromagnetic solution? Best regardsAndrei Postnikov
- ionut ghitiu  a écrit :
>Dear Siesta users,
I've been trying to obtain the band structure of BiFeO3 for some time now, mind 
you it's my first project involving Siesta or dtf so please bear with me. The 
thing is, the band structures that I get are really similar to the ones from 
materialsproject for the same structure (I am using the R3C primitive cell) but 
for the fact that the system seems to be spin degenerate, and so I get no 
splitting of the bands or difference in DOS. 
To give you some details about the system, I'm currently using FR/SR 
pseudopotentials from Pseudo-Dojo and an optimized base, but I also tried with 
the default one (up to TZDP i think). My MeshCutoff usually was 250 Ry, but 
recently changed it to 350 Ry, and I also had a FilterCutoff of 200-250 Ry. 
Also, I've been using a 5x5x5 gamma-centered kgrid and moved to a 7x7x7 . I'm 
initializing the Spin as polarised and setting the spins for the 2 Fe atoms as 
+/- (tried with values and with InitSpin.AF). As tolerances I use 1.d-5 for DM 
and the default one for HTol, but also tried with 1.d-6 for DM and 4.d-4 eV for 
the Hamiltonian. I tried both relaxed and unrelaxed systems. The thing is the 
magnetic moments written in the MullikenPop are close to what I've found in 
literature +-4.37 for Fe, 0.004 for Bi, and 0.029 for O.
Also, I've tried experimenting a little on Fe2O3 with the same R3C structure 
and have the same problem, while on NiO everything went well. For Fe I'm using 
3s2, 3p6, 3d6 and 4s2 electrons for the base.
Could you please give me any hints on what I'm doing wrong or what to try next? 
I feel like I'm going in circles. Thanks a lot!

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Puzzle in Spin Polarized calculation in Siesta

[Andrei Postnikov]

Dear Dmitriy,

as I understand, you are talking about total spin polarisation, which should 
normally
be reproducible throughout the methods. Otherwise, in terms of local spin 
moments,
expressed via Mulliken charges or whatever, the differences in the definition
between methods (or, between basis sets) may substantially add to confusion. 

This said, I'd suggest you to try to separate, as possible sources of worry,
the "local" issues (pseudopotentials, basis) from "global" ones (convergence
in a large system with – probably – narrow bands). Your P.S. gives an 
impression that
you are trying everything possible in different directions at the same time.
For setting the "local" issues, I'd suggest a reasonable compact case, including
your cobalt or whatever, to check against an all–electron method. Here you'd 
check
the magnetism or whatever is essential for your problem to the accuracy needed,
settle the XC etc., and then do not touch these issues anymore.

This settled, the instability of results may come from problematic convergence,
as you have a large system and probably narrow bands which swap back and forth
across the Fermi level. If you have no gap but a metal, then in a large system
you need to be very careful about the dense enough k-mesh, mixing schema 
and weights, electronic temperature etc., to see that your system really 
converges 
and not fluctuates.

Good luck, best regards

Andrei

-- prof. Andrei Postnikov -- 
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
ICPM, 1 Bd Arago - BP 95823, F-57078 Metz Cedex 03, France


- Dmitriy Muzychenko  a écrit :
> Dear Siesta Users and Developers
> 
> I’m really confused by the results of the Spin Polarized calculation in 
> Siesta.
> 
> My system is Bi2Te3 slab with single Co atom substitutions of a Bi atom 
> (typically concentration: 1 Co atom per 490-1100 Bi+Te atoms).
> 
> The confusion arises when I start to analyze the spin properties of this 
> system. Depending on the pseudopotential or/and basis set I observe 
> completely different value of total spin polarization (S_up-S_down). For 
> example, numerous calculations have resulted in the following list of total 
> spin polarization values: 4.33, 4.19, 4.00, 3.99, 3,77, 3.52, 3.00, 2,99, 
> 2,90, 2.63, 2.59, 2.00, 1.99, 1.53 and 0.00. I note, that these values are 
> for exactly the same system, just calculated with a different pseudo(s) or 
> basis set(s).
> 
> Such a wide range of total spin (from 4.33 up to 0.00) raise a logical 
> question, what value can I trust and why? The question is rather not in 
> asking for some “perfect set” of parameters, but in what the general 
> strategy(rule) for analyzing the results should I choose in order to filter 
> out clearly artifact results and then focused on the most reliable one?
> 
> I would be grateful for any help, advise or some “rule” that will allow to 
> filter out incorrect results.
> 
> Thank you in advance for any advice.
> 
> With kind regards,
> Dmitriy Muzychenko
> 
> 
> 
> P.S. I believe that it is useless to attach here any pseudo or basis that I 
> used, since I have tried a lot of possible combinations of pseudo and basis 
> (pseudo and basis: from paper [Comp. Mat. Sci. 98 (2015) 372]; from Siesta 
> Web-Database; my own generated/tested pseudo and basis). I have also allowed 
> Siesta prog to choose basis automatically, as well as used “gen-basis” 
> utility for basis or choose basis manually... Each attempt gives some new 
> value of total spin polarization (perhaps, it can only be noted that the most 
> common numbers are: 4.00, 2.00 and 0.00). Orbital population analysis shows 
> that in all cases polarization is comes mainly from d-shell of Co atom, with 
> the minor extra polarization of surrounding atoms. I have used GGA PBE 
> approximation with Spin Polarization as well as with SOC with the similar 
> results. For the Co species I have tried also semicore (3s2 3p6 3d7 4s2) 
> configuration with close to the same total spin variation from S=4.33 to 0.00 
> depending on basis set or parameters used for pseudopotential generation. 
> Version of Siesta that I used is 4.1.b3/b4 or latest code from “master 
> brunch”.

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Question about the number of electrons considered in the pseudopotential

[Andrei Postnikov]

Dear Tan Shi,
the issue of importance of 3p for Fe has been discussed in 
Eur. Phys. J. D 25, 261–270 (2003) -DOI: 10.1140/epjd/e2003-00209-3
The pseudopotentials I used there were straightforward – generated with r = 
2.00, 2.00, 2.00, 1.50, 0.0, 0.7for both cases, 4 1 0.00and 4 1 6.00
Best regards
Andrei Postnikov
- shi tan  a écrit :
>Dear SIESTA experts,
  For Fe, I noticed that the available pseudopotential from the siesta 
website only includes 4s2 and 3d6 electrons. But in some published works with 
other codes, 3p6 electrons are also included. I suppose that it is more 
accurate with more electrons taken into account. Where can I find those 
pseudopotentials with more electrons or should I construct one by myself? This 
happens to a lot of metals, is it accurate to just use those pseudopotentials 
from the SIESTA website to describe metal behavior? I would really appreciate 
your help! Thank you! 
Sincerely,-- 
Tan SHI

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] PDOS Up and Down

[Andrei Postnikov]

Dear Rayan,in a spin-polarized case, the script writes several columns of PDOS 
(2, or 4 for non-collinear case)into the output file:C --- write down 
accumulated DOS values: -
if (isferm) write (io1,310) efermi
write (io1,304)
do it=1,nt
write (io1,305) ene(it),(dos(it,ispin),ispin=1,nspin)
enddo
close (io1)...Best regardsAndrei
- rayan moukhadder  a écrit :
>Dear siesta users,A successful spin polarized siesta run was done and a .PDOS 
>file has been created. Now I would like to plot PDOS up and down for specific 
>atoms using one of Andrei Postnikov tools which is fmpdos script. 
Using this script I can easily identify n,l,m quantum numbers of the atom but 
how can I specify the z number concerning the spin(up or down), for example if  
I need to plot only the spin up component of the Pz orbital how can I do such 
action.
Thank you.



-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] SCF not converging for Cr2O3

[Andrei Postnikov]

Dear Yuvam:It might be a good idea to fix first the small periodic Cr2O3 with 
less than 720 atoms... then if it works and your cluster don't,I'd say, you 
have problem with your structure (positions of atoms);try to look at them.By 
the way, how about the magnetic moments?An idea to try all different 
pseudopotentials until it starts to convergeis usually not very fruitful...And, 
it might make sense to look not only at the convergence partof the output 
file...Good luckAndrei Postnikov
- Yuvam Bhateja  a écrit :
>Hello everyone,
I was trying to perform CG calculation of Cr2O3 consisting of 720 atoms but it 
is not converging in its first cycle and reaching max scf iterations i.e., 400 
and even 800.I tried lowering the electronic temp to 30 eV but it gets worse.I 
am currently using 50 eV but still not converging. I used both psf and psml 
pseudopotential but still not converging.
I am attaching my input file and a screenshot.Can please anyone help?

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] PDOS at a vacuum on top of a slab

[Andrei Postnikov]

Dear  Soumyajyoti,the easiest approximation to what you probably want 
(especially as you mentionempty spheres) might be to use ghost atoms, see the 
end of the description ofPAO.Basis block in the manual. Notably the atom number 
-100 is supposed to add the basis of spherical Bessel functionsto whenever you 
place your ghost atom, – or, you assign another negative numberand tune the 
basis using all the usual tricks. Then you just look at the PDOSat these ghost 
atoms. (I never tried this, through).Otherwise I'd suggest that you calculate 
LDOS over a sequence of energy intervalsand then integrate each LDOS over 
spacial region you are interested in (a thin layer at some height above the 
surface, in your case).The integrated values, placed along the energy axis, 
will then give you a kind of histogramBest regardsAndrei Postnikov
- Soumyajyoti Haldar  a écrit :
> Dear Siesta Users, > I am trying to calculate PDOS in the vacuum at a Z Å 
> height above a surface. Can you suggest how to proceed with Siesta? In some 
> other codes, I know one can use empty sphere to project DOS. Is there 
> something similar in Siesta one can set up?  If yes can you
> let me know how to do it? > Thanks and regards > Soumyajyoti >Dr. Soumyajyoti 
> Haldar
> Postdoctoral Researcher | AG Heinze
> Institut für Theoretische Physik und Astrophysik
> Christian-Albrechts-Universität zu Kiel> Leibnizstraße 15 | 24098 Kiel | 
> Germany> Email: > hal...@physik.uni-kiel.de


-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] [SIESTA] Performing geometrical optimization

[Andrei Postnikov]

Dear Yuvam, 
if you are new to the code it might be more "robust" to test parts of your 
system (just a single-cell graphene;
separately a single-cell oxide; ... ) before making and relaxing a supercell 
out of them. The point is, the choice of basis functions might be critical in 
Siesta because it represents a deal between accuracy you need and the 
computational effort you would like to spend, and it is easier
to get some feeling about this compromise when testing small systems. 
Exactly for this reason the question "how fast is Siesta", put in general form, 
is ill-defined. The answer would depend
on the accuracy needed (i.e., which properties you are after), moreover – a lot 
– on the geometry of your system,
number of neighbours to each atom, amount of vacuum in the unit cell, etc.
This said, an automatic relaxation of all crystallographic degrees of freedom 
is technically no problem with
MD.VariableCell  set to  T,
but, in a complex system, you might need to fix some of parameters (the lattice 
parameter of the substrate, or else)
before engaging in a global uncontrollable relaxation which might lead to weird 
results. This said, I'd warn you not to expect too much tutorial-like 
introductions from the mailing list(tutorials being accessible in abundance via 
the Siesta web page), but answers to specific questions.Your inquiry seems 
mysterious without indicating who "they" used which scripts for which 
optimization, and what was confusing in it.Best regardsAndrei Postnikov
- Yuvam Bhateja  a écrit :
>Hello everyone,
My name is Yuvam and I am an undergraduate student from Kolkata, India. I am 
new in siesta and have some experience with softwares like Quantum ESPRESSO. I 
want to perform geometrical optimization of my unit cell (atomic position as 
well as cell vectors). I have created my custom made unit cell to accommodate 
the graphene as well as metal oxide for gas sensing. I was following a tutorial 
in which they used scripts for optimization but I find it very confusing and 
very unsuitable as varying all 9 values of lattice vector using loops was very 
cumbersome. Can someone help provide any other method in which system does it 
by its own like in QE? And also, how fast is siesta as compared to other codes 
like QE, as my system consists of 200-300 atoms and using a cluster with 8 
cores and 42 GB RAM.
Thank you in advance.
RegardsYuvam BhatejaB.Tech. 3rd yearE Shibpur

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Restarting a siesta job

[Andrei Postnikov]

- RAHUL SURESH  a écrit :
>Hi..
Thank you for your reply.
My input file is as follows:(...) scf: iscf   Eharris(eV)  E_KS(eV)   
FreeEng(eV)dDmax  Ef(eV)
   scf:1  -557650.096820381.783520381.6763753.10403
   scf:2 -4630269.6330  2077870.1836  2077870.1745864.83778
   scf:3 -2261484.3148   854868.7375   854868.6831*
   scf:4 -1331002.0702   338981.5615   338981.4669425.38740
   scf:5  -719651.9273   106868.1468   106868.1251*
   scf:6  -717147.807168971.176668971.0943426.43790  
Can I just add 
MD.UseSaveCG   .true. 
MD.UseSaveXV  .true.  
DM.UseSaveDM  .true.  

input file without changing anything to the fdf file and re-run?
Thank you
Dear Rahul:you can do whatever you want, but please take 
into account thatyour calculation does not converge and produces strongly 
meaningless results,and this not as it crashed but from the very beginning :   
scf: iscf   Eharris(eV)  E_KS(eV)   FreeEng(eV)dDmax  Ef(eV)
   scf:1   -89709.5924   -89842.8329   -89842.9104  1.93966 -6.9260
timer: Routine,Calls,Time,% =IterSCF 1 152.571  78.20
   scf:2   -89677.0881   -89355.7591   -89355.9226  0.51386 -3.5728
   scf:3   -89796.1277   -89389.5132   -89389.7502  0.97916 -2.3163
   scf:4  -105751.2656   -88983.4222   -88983.7173 54.29636 -8.1884
   scf:5  -491486.7932   -40257.9792   -40258.0558259.16731and in 
fact your output tells you this in plain words:   scf:  999  -779797.2714
78350.029378349.9655683.29828
   scf: 1000  -603514.165749775.908049775.8345*
SCF_NOT_CONV: SCF did not converge in maximum number of steps.
So a good idea might be to fix this, in the first place.Some suggestions:1) Is 
the structure correct? (Did you look at it?)If yes, could a smaller fragment 
than 648 atoms be selectedfor a trial calculation to check that it works 
correctly?2) Did you try a serial calculation? Would it converge? 3) What was 
the mixing parameter?Etc.Please ask this mailing lest responsibly; this is for 
asking questions and discussing specific problems,not for depositing millions 
of lines of output files.Good luckBest regardsAndrei Postnikov

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Restarting a siesta job

[Andrei Postnikov]

Dear Rahul -it depends on where the job was terminated(usually, a more precise 
question leads to a more specific answer). If during the electronic iterations 
- which so far proceeded reasonably - you can normally start from the .DM 
file(see that the switch is set  UseSaveData  T).If during the structure 
relaxation, you can start from the last .XV file.(See that the switch is set  
MD.UseSaveXV  T). It might be a good idea to remove the Broyden history 
files.If during the frozen phonons calculation, you have to startthe 
displacements of the last atom anew, and erase by handthe last added blocks 
(out of 6, calculated for +/- X/Y/Z displacements of each atom)in the .FC file 
-set MD.FCfirst accordingly.Best regardsAndrei Postnikov
- RAHUL SURESH  a écrit :
>Hi users.
The terminal in which the siesta job was running got terminated accidentally. 
How can I restart the job and how far the results are reliable with restart 
data file. 
Thank you
-- 
Regards,Rahul 

-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] Negative translation acoustic modes in phonon band structure

[Andrei Postnikov]

Dear Igor,

your phonons at Gamma are sane (you may check the frequencies,
but otherwise they seem OK). This hints that a priori there is no problem
with the equilibrium geometry.

Now, you run a single cell
> SuperCell_1  0
> SuperCell_2  0
> SuperCell_3  0
so I don't understand how you expect to get phonon dispersions from this.

> (a supercell would otherwise do band folding and band structure would 
> have proportionally more bands).
- right, this is how it is supposed to be. In "SIESTA way" you get force 
constants
from supercell calculation, and then VIBRA Fourier transforms them and
helps you to trace the phonon branches as if in a single cell.
Note that the supercell size is the user's responsibility,
and nothing will automatically prevent you from using a too small one...

If in doubt (and in order to test things), try to make a doubled cell
and calculate "exactly" a zone-boundary phonon...

Best regards

Andrei

-- prof. Andrei Postnikov -- tel. +33-372749149 -- 
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
ICPM, 1 Bd Arago - BP 95823, F-57078 Metz Cedex 03, France




- popov  a écrit :
> Hi,
> 
> My Siesta calculations of phonon band structure of a spinel crystal 
> gives negative-energy
> translational (the lowest two) bands along the Brillouin zone (see the 
> fig. in attachment).
> 
> Suspecting that geometry optimized by CG in Siesta is at a saddle point, 
> before FC run I
> optimized the geometry with various other algorithms, including those 
> from ASE (preconditioned
> LBFGS, mdMin, GPMin) but I could not get rid of negative frequencies for 
> any of the optimization algorithms. The force tolerance for any of these 
> optimizers was 0.005 eV/Ang (default is 0.04).
> 
> The primitive cell is rhombohedral and I use vibra for postprocessing of 
> the FC calculations.
> Input for Vibra is similar to fdf for FC calculations, with addition of:
> 
> SuperCell_1  0
> SuperCell_2  0
> SuperCell_3  0
> AtomicDispl  0.06  Ang
> Eigenvectors True
> 
> and I added masses of atoms in %block AtomicCoordinatesAndAtomicSpecies 
> in 5th column
> and specified AtomicCoordinatesFormat  NotScaledCartesianAng.
> (should I consider InitMesh instead?)
> corresponding to the real space grid step division of 5.9114/130 = 0.045 
> Ang, where 5.9114 Ang
> is the module of lattice vectors. I used MD.FCDispl 0.06 Ang.
> 
> I use a primitive cell with 14 atoms for the FC calculations, i.e. I did 
> not make a supercell
> in order to reconstruct 14*3=42 bands in band structure that is common 
> in literature for spinels
> (a supercell would otherwise do band folding and band structure would 
> have proportionally more bands).
> 
> Now I am running some additional calculations with finer MeshCutoff and 
> MD.FCDispl (1400 Ry, 0.02 Ang).
> 
> I am out of options. Do you have a clue what could be cause of imaginary 
> frequencies of the
> translational modes throughout the BZ and how to get rid of them?
> 
> Thanks in advance,
> Igor
> 
> 


-- 
SIESTA is supported by the Spanish Research Agency (AEI) and by the European 
H2020 MaX Centre of Excellence (http://www.max-centre.eu/)

Re: [SIESTA-L] problem

[Andrei Postnikov]

Dear doaa haleem :
in your “111” case, your unit cell is twice longer in the X direction than in 
other two(did you have a look at the structure figure??)
In other two cases, the lattice vectors are more cubic, yet you have Hf in the 
“A” positionand Ba in the “B” position (in the oxygen octahedron) of the 
perovskite. Is this as supposed to be?Won’t this be the other way around?
Concerning the spin different from 0, this might be due to settings in your 
LDAU block.You switch some U on Hf and set SpinPolarized =T; what do you 
expect, in fact?
In order to better see what’s going on you can add some printout (Mulliken 
charges…).I’d say, try first to get a reasonable structure in non-magnetic case 
without LDA+U...
Best regards
Andrei Postnikov
- doaa haleem  a écrit :
>I am writing this mail to discuss a problem I have faced when I was 
>calculating using siesta package, I am using siesta 4.0 , I am working on  
>some peroveskite structure , at the beginning I have used calculations using 
>GGA only ,it is expected that the final structure (XSF) to be  the same as the 
>start structure (FCC) ,but it isn't and it consists only from 5 atoms with 
>wrong positions..(I will attach the fdf file, XSF, and gnuplot draw that 
>resulted for the previous case named after "case 1") The same have happened 
>when using GGA+U (case 2)Trying to solve this problem I have used the geometry 
>constraint to fix Hf atom and all  3-O atoms as in the perovskite BaHfO3, the 
>resulted structure is wright but -in the DOS and PDOS - the spin up DOS isn't 
>the same as the spin down DOS ( I will send the fdf file ,  XSF and   gnuplot 
>draw for the previous case named after "case ") I tried also to use "variable 
>cell command , repeating the run and changing the k-grid and mesh cut off but 
>the problem is still existing .It is notable that I am using the pseudo 
>potentials uploaded on siesta web page So, please help me to solve this 
>problem .I need the the structure to be right, the spin up the same as the 
>spin down and the total spin is zero .Thank you 

Re: [SIESTA-L] EXTRACTION OF PDOS --> modified version of fmpdos

[Andrei Postnikov]

Dear Siesta users,in view of some modifications of the PDOS format in recent 
versions of Siesta,my script fmpdos.f (for extraction of atom-resolved partial 
densities of states)required modification. The updated version (hopefully able 
to manage bothold and new formats) is in attachment, and also available 
athttp://www.home.uni-osnabrueck.de/apostnik/Software/fmpdos.fThank you very 
much in advance for testing and sending me bug reports. Best regardsAndrei 
Postnikov- James Sifuna  a écrit :
>Hello, 
I have successfully done the dos and pdos calculations, 
I need to extract the pdos of atom 1 with n=4, l=0, m=0 etc. but i get an 
error..

siesta-4.0.2/Util/Contrib/APostnikov/fmpdos 
  Input file name (PDOS):
SrTiO3.PDOS
  Output file name :
Sr.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):
1
  Extract data for n= ... (0 for all n ):
0
  Unknown identifier in PDOS file, line6
  -0.69252E+02  
 


Best wishes,
J.
	program fmPDOS
!
!   extracts partial densities of states from PDOS file of SIESTA,
!   taking care of flexible format in index, atom_index, etc. fields
!   Andrei Postnikov, University Paul Verlaine - Metz (Jul 2010)
! University of Lorraine  (Jan 2019)
!   andrei.postni...@univ-lorraine.fr
!
  implicit none
  integer ii1,io1
  parameter (ii1=11,io1=12)
  integer nt,nmax,i0,ispin,nspin,norbs,it,npts,nene,is,idos,l,lref,
 .m,mode,mref,zeta,zvalue,atind,n,nref,indref,index,nline,
 .iquoted,iparsed
  parameter (nmax=1)   ! max. number of energy values in DOS
  double precision ene(nmax),dos(nmax,4),dos1(4),dparsed,efermi
  character inpfil*60,outfil*60,string*80,llabel*80,rlabel*80,
 .  species*6,squoted*6,chlab*6,owrite*1,s1*1,polar*6,
 .  dtail*7
  logical filexist,addos,nptdef,ene_fo,isferm
  external iquoted,squoted,iparsed,dparsed

  do it=1,nmax
  do is=1,2
dos(it,is)=0.d0
  enddo 
  enddo 
  nptdef = .false.   ! We don't know whether the file contains 
  isferm = .false.   ! Show the Fermi energy if included in the PDOS

! --- open input and outpute files: --
  write (6,*) ' Input file name (PDOS):'
  read (5,*) inpfil
  inquire (file=trim(inpfil), exist=filexist)
  if (filexist) then
open (ii1,file=trim(inpfil),status='old',form='formatted')
  else
write (6,*) ' File does not exist!'
stop
  endif

  write (6,*) ' Output file name :'
  read (5,*) outfil
  inquire (file=trim(outfil), exist=filexist)
  if (filexist) then
write (6,*) ' File ',trim(outfil),' exists. Overwrite? (Y/N)'
read (5,*) owrite
if (owrite.eq.'Y'.or.owrite.eq.'y') then
  open (io1,file=trim(outfil),form='formatted',status='REPLACE')
else
  write (6,*) ' Then rename is first. Bye...'
  stop
endif
  else
open (io1,file=trim(outfil),form='formatted',status='NEW')
  endif
  write (io1,309) trim(inpfil)

! --- select projected DOS which are needed: -
  100 continue
  write (6,*) ' Extract data for atom index',
 .' (enter atom NUMBER, or 0 to select all),' 
  write (6,*) ' or for all atoms of given species',
 .' (enter its chemical LABEL):'
  read (5,*) string
! check first character of string, whether it is decimal
  i0 = ichar(string(1:1))
  if (i0.ge.48.and.i0.le.57) then
!   1st character passed is a number;
!   assume the whole string is atom number
mode = 1
read (string,*,err=401) indref
write (io1,301)indref
  elseif ((i0.ge.65.and.i0.le.90).or.(i0.ge.97.and.i0.le.122)) then
!   1st character passed is a character [A-Z,a-z];
!   assume the whole string is atom label
mode = 2
read (string,*,err=402) species
write (io1,302) trim(species)
  else
!   1st character is a special symbol; presumably an error
write (6,306) i0
goto 100
  endif

  write (6,*) ' Extract data for n= ... (0 for all n ):'
  read (5,*) nref
  if (nref.ne.0) then
write (6,*) ' Extract data for l= ... (-1 for all l ):'
read (5,*) lref
if (lref.ne.-1) then 
  write (6,*) ' Extract data for m= ... (9 for all m ):'
  read (5,*) mref
else
  mref=9
endif 
  else 
lref=-1
mref=9
  endif 
  write (io1,308) nref,lref,mref

  ene_fo = .false.  ! are we in the energy reading mode ?
  addos  = .false.  ! add this DOS to the result ?
  nene=0! counter for energy values
  nline = 0 ! global conter for lines in PDOS file

! --- read and analyze the PDOS file line by l

Re: [SIESTA-L] DOS calculation with fermi shifted to zero

[Andrei Postnikov]

Dear Riya,the DOS / PDOS come out in the "native" energy scale of the code(in 
which Fermi energy - listed in the last column of electronic iteration 
printout, under siesta: iscf   Eharris(eV)  E_KS(eV)   FreeEng(eV)   dDmax  
Ef(eV)- is typically several eV in the negative).In order to plot, you shift it 
"by hand" within your favourite plotting program.Thank you for your attention 
to Javier Junquera's tutorial; there is no special "Andrei Postnikov method" to 
plot the density of states...
Best regardsAndrei Postnikov
- Riya Rogers  a écrit :
>Dear All,
I am calculating DOS of my system using the Andrei Postnikov method (given in 
tutorial below) :
http://personales.unican.es/junqueraj/JavierJunquera_files/Metodos/Structuralproperties/PDOS-SrTiO3/Exercise-PDOS-SrTiO3.pdf

I have a doubt that whether this method gives the DOS with fermi automatically 
shifted to zero?
or 
Do we have to manually shift fermi to zero?
Please tell me how to get DOS with fermi shifted to zero, I will be highly 
obliged.
Thanking you
RegardsRiya

Re: [SIESTA-L] File Not Converging

[Andrei Postnikov]

Dear Harsh Shah,it may sound banal, but - if you hope for a useful answer you 
should outline your problem in a meaningful way. "A certain system with 110 
atoms did converge but the other one with 769does not" - this is not much of 
information.You define 2000 SCF iterations and 5000 CG steps - who of them do 
not converge?What is your mixing parameter, anyway? How does your failed 
convergence look like?
In some systems (no idea what is yours) with many nearly degenerate statesa 
very drastic reduction of mixing might help(to get the convergence trend right, 
before trying more efficient mixing schemes).Best regards
Andrei Postnikov
- Harsh Shah  a écrit :
>Dear Siesta User, 

I am running a file in parallel with around 769 atoms in stampede2 
supercomputer and my file is not converging.
It has been running without giving any error for around 10-12 days and I think 
it should not take that much time.I ran few other files having 110 atoms and it 
is giving proper output.
I am using:MD.MaxForceTol0.025MeshCutoff(in Ry)  80MaxSCFIterations   
2000MD.NumCGSteps   5000DM.Tolerance 1.d-4
Is there anything which I should change from above?
What can be other possibilities do to which the file is not converging?
I would really appreciate if you could help me out.
Regards,Harsh Shah

Re: [SIESTA-L] atom resolved phonon density of states

[Andrei Postnikov]

Dear Karolina:thanks for your interest; please try the attached code(a single 
fortran source file, including all the subroutines called).Just compile it.It 
is not separately documented, but you just run it and answer the 
questions.Please ask me if something is not clear or goes wrong,best 
regardsAndrei-- prof. Andrei Postnikov -- tel. +33-372749149 -- University of 
Lorraine - Laboratoire de Chimie/Physique - A2MCICPM, 1 Bd Arago - BP 95823, 
F-57078 Metz Cedex 03, 
France
- Karolina Milowska  a écrit :
>Dear All,

I would like to plot a phonon DOS resolved over particular atoms in my system, 
in a similar fashion to the standard projected DOS. How may I extract this 
information in practice?

I have found  a presentation from 2010:
http://www.home.uni-osnabrueck.de/apostnik/Lectures/APostnikov-Phonons.pdf
which shows what I would like to get. On slide 21, there is an information 
about utility - phdos. Where can I find it? 


Kind regards,
Karolina

C
C   phdos, a script to calculte phonon density of states
C  in a supercell,
C  using phonon eigenvectors provided by "vibrator"

C  Written by Andrei Postnikov, Oct 2005   Vers_0.2
CXV read format bug fixed   May 2006
C  apost...@uos.de
C
  program phdos
  implicit none
  integer ii1,ii2,io1,io2,ivmin,ivmax
  parameter (ii1=11,ii2=12,io1=14,io2=15)
  integer ialloc,iat,nat,iatmin,iatmax,mode,nodif,idif,npoints
  integer, allocatable :: ityp(:),iz(:),jat(:,:),jtyp(:),nna(:)
  double precision tau(3,3),qvec(3),qq(3),delta
  double precision, allocatable :: mass(:),coor(:,:),
 .   freq(:),evr(:,:,:),evi(:,:,:),wwsum(:,:),zz(:)
  character inpfil*60,outfil*60,syslab*30,qlab*1
  character*2, allocatable :: label(:)
  external test_xv,read_xv,read_ev,full_dos,qres_dos
C
C string manipulation functions in Fortran used below:
C len_trim(string): returns the length of string
C   without trailing blank characters,

C --- read lattice vectors and atom positions from the .XV file:
  write (6,701)
  701 format(' Specify SystemLabel of .XV file: ',$)
  read (5,*) syslab
  inpfil = syslab(1:len_trim(syslab))//'.XV'
  open (ii1,file=inpfil,form='formatted',status='old',err=801)
  call test_xv(ii1,nat)
  allocate (ityp(nat))
  allocate (iz(nat))
  allocate (mass(nat))
  allocate (label(nat))
  allocate (jtyp(nat))
  allocate (nna(nat))
  allocate (jat(nat,nat))
  allocate (coor(1:3,1:nat),STAT=ialloc) 
  if (ialloc.ne.0) then
write (6,*) ' Fails to allocate space for ',nat,' atoms.'
stop
  endif
  call read_xv(ii1,nat,ityp,iz,tau,mass,label,coor)
  close (ii1) 
C --- count different types in the XV file. Note that they can be changed
C arbitrarily (by hand) in the XV file  as compared to the Siesta run. 
C E.g. to select some atoms of type=4 as different, assign them type=14.
C They will then get a separate column in output files.
  nodif=0
  do 98 iat=1,nat
if (nodif.gt.0) then
  do idif=1,nodif
if (ityp(iat).eq.jtyp(idif)) then  !  add atom to existing type
   nna(idif)=nna(idif)+1
   jat(idif,nna(idif))=iat
   goto 98
endif
  enddo
endif
nodif=nodif+1 !  add a new type
jtyp(nodif)=ityp(iat)
nna(nodif)=1
jat(nodif,1)=iat
   98 continue
  write (6,"(i5,' different types found in the XV file:')") nodif
  do idif=1,nodif
write (6,"(i5,' atoms of type ',i5)") nna(idif),jtyp(idif)
  enddo
C --- read and store frequencies /eigenvectors from the .vector file,
C q =(0 0 0)  only :
C allocate for all modes, 1 through nat*3 :
  ivmin=1
  ivmax=nat*3
C (change ivmin, ivmax above if want to restrict to less modes)
  allocate (freq(ivmin:ivmax))
  allocate (evr(1:3,1:nat,ivmin:ivmax))
  allocate (evi(1:3,1:nat,ivmin:ivmax))
  allocate (zz(1:nodif))
  allocate (wwsum(1:nodif,ivmin:ivmax),STAT=ialloc)
  if (ialloc.ne.0) then
write (6,*) ' Fails to allocate space for vibration modes'
stop
  endif
  write (6,702)
  702 format(' Specify SystemLabel of .vectors file: ',$)
  read (5,*) syslab
  inpfil = syslab(1:len_trim(syslab))//'.vectors'
  open (ii2,file=inpfil,form='formatted',status='old',err=801)
  call read_ev(ii2,nat,qq,ivmin,ivmax,evr,evi,freq)
  write (6,*) 'opened and read ',inpfil
  close (ii2)
C ---
  write (6,*) 'You have two options: Total density of states ',
 .'(sum of squares of eigenvectors),'
  write (6,*) 'or Q-projected density of states ',
 .'(convolution wit

Re: [SIESTA-L] VACF and PHONON DOS in SIESTA

[Andrei Postnikov]

Dear Suman:
you can tryhttp://www.home.uni-osnabrueck.de/apostnik/download.html  -> 
VeloCorrF
Thank you for reporting bugs,   Sincerely
Andrei Postnikov
- Suman Chowdhury  a écrit :
>Dear SIESTA user : Is there any utility program is SIESTA which can calculate 
>the velocity autocorrelation function and the fourier transform of it, because 
>the fourier transform of velocity autocorrelation function gives the phonon 
>density of states which I want to calculate..

Re: [SIESTA-L] Error in running grid2cube

[Andrei Postnikov]

- Ankita Joshi  a écrit :>Dear Siesta Users,
While plotting LDOS for my system using a post-processing tool grid2cube, I 
encountered following error:STOP grid2d: Parameter MAXP too small
Does anyone knows how to solve this problem? Please give your suggestions. 


Dear Ankita:This is not an error but a security message to be taken 
literally.The suggestion to solve this problem is,open  grid2cube.fin your 
favourite editor, correct the line   parameter (maxp = 1200)
to a larger value and recompile. You could have found the solution 
yourselfbrowsing the grid2cube.f  for the stop message:
  np = mesh(1) * mesh(2) * mesh(3)  if (np .gt. maxp) stop 'grid2d: 
Parameter MAXP too small'
( Do you really need the mesh more dense than that? )Best regards
Andrei Postnikov

Re: [SIESTA-L] cell vector not converged

[Andrei Postnikov]

- Najmeh Honari  a écrit :
Dear Andrei,
In terms of converging our simulations for k-mesh values, how do we know
> whether the fluctuations of the "total energy" and 
> 
> "Cell Vector Mod"
> are big or small? 
Dear Najmeh,do you refer to the recently posted results of Sunetra Das (are 
they "your simulations"?),or to something else? - Just in order to define what 
we are talking about.I'd say, the general approach is like this:- you converge 
the total energy (differences, see my previous comment) vs the k-ponts to the 
accuracy you need for your particular problem. - once you decide that the 
convergence is good, you write down the lattice parameters.To look at how do 
the lattice parameters fluctuate on increasing the k-mesh densitydoes not make 
much sense, because the lattice parameter is not a variational property.  Best 
regardsAndrei Postnikov

On 14 July 2018 at 03:02, Andrei Postnikov  
wrote:
Dear Sunetra,
1) your system is 2D; your "Cell Vector Mod" is presumably the lattice 
parametouler in (x,y);how about the third parameter? Is there a chance that you 
have not fix it so that itfluctuates under the effect of stress tensor 
component which is a numerical noise?
2) "Is the free energy convergence enough?" - the answer fully depends onyour 
objective; normally this is not the absolute value of energy that you should 
care about -it has no special meaning anyway - but the energy difference 
between cases(geometries, magnetic solutions, etc.) which you want to compare 
in your study.
Best regards
Andrei Postnikov


- Sunetra Das  a écrit :
>Dear experts,
I am trying to converge my system unitcell structure for k - mesh values. Even 
though the Free Energy value has converged for less denser mesh, the unit cell 
parameter has not converged for even 21x21x1 k mesh. The system is a 2D one 
with 3 atoms in the unit cell. I am attaching the file containing the k values 
and the corresponding free energy and cell parameter values along with the 
plotted graph. 

please have a look and help me. Is the free energy convergence enough? I have 
used systems before where I have checked if both the parameters have converged 
for a given k mesh.
Thank you.
Regards,Sunetra Das

Re: [SIESTA-L] cell vector not converged

[Andrei Postnikov]

Dear Sunetra,
1) your system is 2D; your "Cell Vector Mod" is presumably the lattice 
parametouler in (x,y);how about the third parameter? Is there a chance that you 
have not fix it so that itfluctuates under the effect of stress tensor 
component which is a numerical noise?
2) "Is the free energy convergence enough?" - the answer fully depends onyour 
objective; normally this is not the absolute value of energy that you should 
care about -it has no special meaning anyway - but the energy difference 
between cases(geometries, magnetic solutions, etc.) which you want to compare 
in your study.
Best regards
Andrei Postnikov
- Sunetra Das  a écrit :
>Dear experts,
I am trying to converge my system unitcell structure for k - mesh values. Even 
though the Free Energy value has converged for less denser mesh, the unit cell 
parameter has not converged for even 21x21x1 k mesh. The system is a 2D one 
with 3 atoms in the unit cell. I am attaching the file containing the k values 
and the corresponding free energy and cell parameter values along with the 
plotted graph. 

please have a look and help me. Is the free energy convergence enough? I have 
used systems before where I have checked if both the parameters have converged 
for a given k mesh.
Thank you.
Regards,Sunetra Das

Re: [SIESTA-L] Problem in SCF convergence in spin polarized calculations

[Andrei Postnikov]

Dear Ananya:OK it doesn't converge, but are the intermediate (not converged) 
results meaningful, or crazy?Look at the Mulliken charges - are they more or 
less stable (and reasonable),in the course of iterations, or wildly 
fluctuating?Where is the magnetism - where it belongs (Fe impurity? ribbon 
edges?) Did you initialised the magnetism correctly?Do you see in the density 
of states more or less what you expect(knowing the DOS of pure graphene and 
knowing how the Fe 3d statesare supposed to be placed)?Etc. etc.So may things 
to look at, beyond just the convergence...My wild guess would be, ALL your 
carbon atomsare initialised to maximal spin, and then the things go astray...
Best regardsAndrei Postnikov
- Ananya Rajpoot  a écrit :
>
> 
>Dear users,I am interested in simulation of Graphene ribbons doped with Fe. 
>When I proceed with the spin-polarized calculations, the SCF in first CG run 
>didn't converge. It is taking 500 SCF cycles or even more during 
>convergence.Should I change pulay parameters? Or Linear mixing schemes? If 
>pulay then should I increase weight or number of pulay??Also, what else should 
>I correct to get faster convergence?? or Electronic temp should be reduced?
My fdf file parameters are as follows:PAO.BasisSize DZPXC.functional 
GGAXC.authorsPBE%block kgrid_Monkhorst_Pack
   1   0   0  0.0   0   1   0  0.0   0   0  49  0.0%endblock 
Kgrid_Monkhorst_PackMeshCutoff 370.0 Ry
MaxSCFIterations 5000
SCF.Mixer.Method PulaySCF.Mixer.Variant originalSCF.Mixer.Weight 
0.1SCF.Mixer.History 2SCF.Tolerance.DM 1.d-4
MD.TypeOfRun CGMD.MaxForceTol = 0.01 eV/AngMD.NumCGsteps 5000MD.VariableCell 
.false.
SCF.Mix.First falseSCF.Mix.AfterConvergence falseSCF.RecomputeHAfterScf 
falseSCF.Mix densitySCF.Converge.DM true
Spin polarized
SCF.Mix.Spin allFixSpin falseDM.InitSpinAF false
WriteCoorXmol trueSolutionMethod   diagon
Thank u so much...With Regards
Anaya
> 
>

Re: [SIESTA-L] Changing the crystal direction of Si

[Andrei Postnikov]

Well, your unit cell is now 4 times larger,so you need to specify 8 atoms per 
cell...If this is really what you want...Whatever is the meaning of your 
"change the crystal direction to <100> "...(do you want to expand the supercell 
along [100] for tracing the phonon dispersion?A doubled supercell [ 1 0 0 ] [ 0 
1/2 -1/2 ] [ 0 1/2 1/2] with 4 Si atoms should work fine)Best regardsAndrei 
Postnikov
- Alaa Akkoush (Student)  a écrit :
>> I have build the Silicon which has diamond structure and calculate it's 
>> dispersion relation.
> I set the lattice vectors those of fcc and place two silicon atoms one at 0 
> and other at (0.25,0.25,0.25)
> as follows:

> %block LatticeVectors
> 0.0 0.5 0.5
> 0.5 0.0 0.5
> 0.5 0.5 0.0
> %endblock LatticeVectors

> %block AtomicCoordinatesAndAtomicSpecies
>  0.000.000.00   1   28.0855
>  0.250.25   0.251   28.0855
> %endblock AtomicCoordinatesAndAtomicSpecies
> I want to CHANGE THE CRYSTAL DIRECTION TO <1 0 0>, so i changed the lattice 
> vector to 

> %block LatticeVectors
> 1.0 0.0 0.0
> 0.0 1.0 0.0
> 0.0 0.0 1.0
> %endblock LatticeVectors
> But the system is unstable!!!> Any help? 

Re: [SIESTA-L] Fwd: problem with converging mesh-cutoff and K-grid points

[Andrei Postnikov]

Dear Sunetra,please have a look 
intohttp://www.home.uni-osnabrueck.de/apostnik/Lectures/SIESTA-tuto.pdfin 
particular from page 10 on.Best regardsAndrei Postnikov

- Sunetra Das <sunetra.das...@gmail.com> a écrit :
>Dear Andrei Postnikov Sir,
Can you please guide me in how to do the 'eggbox test' ? I have gone through 
the tutorials and examples but am still unclear on how to do that test for a 
system. Will be very much thankful if you let me know the steps on testing the 
eggbox effect and eliminating it. I have seen the term 'eggbox energy' in 
siesta output file after completion of a run, but confused what to decipher 
from it. Thanking you and awaiting a reply. 
Warm regards and wishing you a happy new year in advance,Sunetra Das.

Re: [SIESTA-L] Fwd: problem with converging mesh-cutoff and K-grid points

[Andrei Postnikov]

Dear Najmeh,you are worrying about z-components of forcesand try to help it by 
increasing the k(xy) mesh. It means that your bands dispersion becomes ever 
better in the plane,however, I don't think this would have an important effect 
on the z-force! 
Three observations:1) Do you really have a problem?Your z-forces on atoms 1 and 
2 do nicely sum up to zero within 0.001 - that is fine, -
and the fact that they (individually) fluctuate with k simply reveals the fact 
that your structure is not at equilibrium. Try to get a better equilibrium - 
see 2).If you ultimate goal is phonons, why don't you calculate just Gammaand 
check if you get the acoustic modes close enough to zero?2)In order to suppress 
fluctuations of forces / energy, 
the MeshCutoff is probably much more important that k-mesh.Apparently you tried 
300 Ry and did not increase it any further.It might be not sufficient, 
depending on pseudos / basis.(However, it might well be OK in your case, see 
1). In order to understand what's going on, the eggbox test (rigidly moving all 
the atoms... along Z, your "problematic" direction) for different values of 
MeshCutoff would help.With insufficient MeshCutoff, the increase of the k-mesh 
alone won't help you.(and... do you really need 20 Ang in the z-direction? You 
only needso much that the basis functions won't overlap).3)The k-mesh is not a 
variational parameter,
so you cannot expect a systematic lowering of the total energy(although it 
usually happens, albeit with fluctuations).For the phonon calculations, the 
total energy as such (i.e.,its absolute value) is irrelevant anyway.
Best regardsAndrei Postnikov

- Najmeh Honari  a écrit :
>Dear
> Siesta users,

> 
>I’m
> trying to simulate 2D borophene δ6, in order to calculate phonon dispersion 
> of the material. I did convergence
> test for mesh-cut-off and k-points. I selected mesh cutoff=300Ry and started
> k-grid from 12x12x1 and increased it to 40x40x1. The issue is that neither 
> “total
> energy” nor the “forces on atoms” are converged in other words fluctuations of
> energies and forces are obsereved. Here are a table of some energy and forces
> for  different k-grid meshes of borophene
> structure:
> 
>>
>K-grid
>
>Total
> Energy (eV)
>
>Atomic
> forces (eV/A)
>
>15x15
>
>-206.414196
>
>1   -0.01   -0.05  
> -0.131499
>2   -0.02   -0.040.131227
>
>17x17
>
>-206.410794
>
>10.00   -0.03  
> -0.423281
>20.04   -0.010.423056
>
>20x20
>
>-206.413284
>
>10.01   -0.06  
> -0.287325
>2   -0.020.020.287087
>
>23x23
>
>-206.414429
>
>10.00   -0.00  
> -0.355986
>2   -0.02   -0.050.355733
>
>25x25
>
>-206.423565
>
>1   -0.01   -0.00  
> -0.416997
>2   -0.02   -0.030.416772
>
>27x27
>
>-206.424037
>
>1   -0.02   -0.06  
> -0.308230
>2   -0.010.000.307990
>
>30x30
>
>-206.422290
>
>1   -0.01   -0.04  
> -0.313348
>2   -0.00   -0.020.313111
>
>33x33
>
>-206.423341
>
>1   -0.02   -0.08  
> -0.368439
>2   -0.01   -0.050.368198
>
>35x35
>
>-206.422301
>
>10.02   -0.04 
>  -0.275519
>2   -0.020.010.275263
>
>38x38
>
>-206.420742
>
>10.020.00  
> -0.288896
>2   -0.000.010.288655
>
>40x40
>
>-206.420015
>
>1   -0.01   -0.05  
> -0.357026
>2   -0.000.000.356782
>
> 
> 
> 
>The
> input.fdf file that is used in simulations is as follow:
> 
> 
> 
>
>SystemName  B
>SystemLabel Borophene
>NumberOfAtoms   2
>NumberOfSpecies 1
>##
>#Basis#
>##
>PAO.BasisSize   DZP  
>##
>#ChemicalSpecies#
>##
>%block
> ChemicalSpeciesLabel
>  1  
> 5   B
>%endblock
> ChemicalSpeciesLabel
>##
>LatticeConstant   1.00 Ang
> 
>%block LatticeVectors
>1.6   0. 
> 0.
>0.0   2.8200 
> 0.
>0.0   0. 
> 20.000
>%endblock
> LatticeVectors
> 
>AtomicCoordinatesFormat
> Ang
>

>%block
> AtomicCoordinatesAndAtomicSpecies
>0.0 0.0 0.0
> 1
>0.8 1.41000 
>   0.91000 1
>%endblock AtomicCoordinatesAndAtomicSpecies
>

>BandLinesScale  
>ReciprocalLatticeVectors
>

>%block BandLines   
>1 0.00 0.00 0.00 Gamma
>200   0.50 0.00 0.00 X
>200   0.50 0.50 0.00 S
>200   0.00 0.50 0.00 Y 
>200   0.00 0.00 

Re: [SIESTA-L] Structural relaxation of 2D carbon atom chain

[Andrei Postnikov]

Dear Tomasz, your carbon chain presumably doesn't want to bendbecause you keep 
the lattice constant (along the chain) too large, and fixed.Otherwise, if you 
want help please specify your problem - what you tried and what goes wrong.<< 
Problem starts when after hundred of calalculations the structer remaind the 
same >>does not seem to give a clue.Otherwise,you may have a look into my 
tutorial on a similar issue(an atomic - gold - chain; tendency to bend; 
stability; phonons) -an explaining presentation and the work directory with 
guidelines,http://www.home.uni-osnabrueck.de/apostnik/Lectures/Exercis-Phonons.pdfhttp://www.home.uni-osnabrueck.de/apostnik/Lectures/Execs_Barcelona_2017.zip
at a Siesta school in Barcelona, May 2017.Good luckAndrei- 
-- prof. Andrei Postnikov -- tel. +33-372749149 
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
ICPM, 1 Bd Arago - BP 95823, F-57078 Metz Cedex 03, France



- Tomasz Kepczynski <kepczynskisie...@gmail.com> a écrit :
>Ive just started using Siesta and wanted to simulate something basic. My tutor 
>adviced me to do for example carbon atomic chain but asked me to try simulate 
>it with 2 carbon atoms in a primitive cell with one placed out of the Main 
>axis So the atomic chain would be bend/ VVV shaped. Problem starts when after 
>hundred of calalculations the structer remaind the same and creats carbon 
>"dimers" at best Hope you can help me start my adventure with Siesta

Re: [SIESTA-L] Doping

[Andrei Postnikov]

- masoudeh maleki <masodeh.mal...@gmail.com> a écrit :
>Dear Andrei,Thank you so much for your email. I'm so confused because I'm 
>trying to dope tinoxide with nitrogen by siesta code and I can’t optimize 
>three parameters such as k point, mesh cut off and lattice constant. I mean 
>that plot of the total energy against these parameters isn’t as well as pure 
>tinoxide. Before doping I optimize these parameters in bulk tinoxide and then 
>relaxed and after relaxed the structure. Should I do all these process in 
>doped tinoxide again? And would you please kindly guide me what should I do to 
>have acceptable plots in doped case? Should I change the coordinates doped 
>nitrogen atom manually or subsituating nitrogen instead of oxigen or tin is 
>enough ?
Yours SincerelyMasoudeh Maleki 
Dear Masoudeh -indeed, it looks like a big confusion. OK, slowly: > I can’t 
optimize three parameters such as k point, mesh cut off and lattice constant.
These "parameters" are not in the same standing. The lattice constant is your 
result,you can indeed optimise it by varying it and looking for total energy 
minimum. The mesh cutoff and k-points a responsible for the accuracy of 
calculation.There are no "optimal" value of them that would minimise the total 
energy.Typically, as you increase them the total energy will consistently go 
down.The larger (mesh cutoff and k-points cutoff) the better; however you won't 
like to overkill.Therefore you make a reasonable free choice of these 
parameters driven by analysis of whatis the accuracy of the total energy you 
really need. This is never the total energyby itself, but a difference of 
energies between different states of phases.As the target accuracy may 
considerably vary from task to task, so may varythe technical parameters of 
your calculation.
> Should I do all these process in doped tinoxide again? If the chosen mesh 
> cutoff value is OK for pure SnO2-and- also for some system containing the 
> dopant it ought to be OKfor the SnO2 + dopant as well.The "necessary" k-mesh 
> essentially depends on the supercell size.Keep about the same density of 
> k-points (i.e. the same real spacecutoff radius) as you increase your 
> supercell. 
> Should I change the coordinates doped nitrogen atom manually > or 
> subsituating nitrogen instead of oxigen or tin is enough ?Without knowing 
> what is your task nobody will tell you what is "enough".In principle, the 
> lattice will relax around a substitutional impurity.What's a big problem to 
> let it , in a Siesta calculation?
Best regards
Andrei Postnikov

On Wed, Sep 27, 2017 at 11:31 PM Andrei Postnikov 
<andrei.postni...@univ-lorraine.fr> wrote:


- masoudeh maleki <masodeh.mal...@gmail.com> a écrit :
>Dear sir,
I'm using siesta code to investigate the effect of doping on bandstructure and 
PDOS of periodic crystal of tin oxide with tetragonal unitcell. I made a 2*2*2 
supercell of it, and replaced one O atom with my dopant. According to the 
siesta manual, can l use NetCharge and simulatedoping together for my bulk 
structure to simulate doping? Or it is impossible? And if it is incorrect how 
can I dope my bulk structure in siesta?
Yours sincerelyMasoudeh maleki
Dear Masoudeh,
if you explicitly add the dopant atom to your system, you don't need any 
NetChargesince your system is not charged. (The dopant atom comes with its 
correct nuclear chargeand number of electrons). As explained in the manual, the 
NetCharge option is providedto treat charged finite systems (molecules). Best 
regards
Andrei Postnikov

Re: [SIESTA-L] Doping

[Andrei Postnikov]

Dear Masoudeh -indeed, it looks like a big confusion. A pity because these 
issues were so many timesaddressed in tutorials and - I am sure - in the 
mailing list.OK, slowly: > I can’t optimize three parameters such as k point, 
mesh cut off and lattice constant.
These "parameters" are not in the same standing. The lattice constant is your 
result,you can indeed optimise it by varying it and looking for total energy 
minimum. The mesh cutoff and k-points a responsible for the accuracy of 
calculation.There are no "optimal" value of them that would minimise the total 
energy.Typically, as you increase them the total energy will consistently go 
down.You just make a free choice of these parameters driven by analysis of 
whatis the accuracy of the total energy you really need.
> Should I do all these process in doped tinoxide again? If the given mesh 
> cutoff value is OK for pure SnO2 -and- also for

- masoudeh maleki <masodeh.mal...@gmail.com> a écrit :
>Dear Andrei,Thank you so much for your email. I'm so confused because I'm 
>trying to dope tinoxide with nitrogen by siesta code and I can’t optimize 
>three parameters such as k point, mesh cut off and lattice constant. I mean 
>that plot of the total energy against these parameters isn’t as well as pure 
>tinoxide. Before doping I optimize these parameters in bulk tinoxide and then 
>relaxed and after relaxed the structure. Should I do all these process in 
>doped tinoxide again? And would you please kindly guide me what should I do to 
>have acceptable plots in doped case? Should I change the coordinates doped 
>nitrogen atom manually or subsituating nitrogen instead of oxigen or tin is 
>enough ?
Yours SincerelyMasoudeh Maleki 
On Wed, Sep 27, 2017 at 11:31 PM Andrei Postnikov 
<andrei.postni...@univ-lorraine.fr> wrote:


- masoudeh maleki <masodeh.mal...@gmail.com> a écrit :
>Dear sir,
I'm using siesta code to investigate the effect of doping on bandstructure and 
PDOS of periodic crystal of tin oxide with tetragonal unitcell. I made a 2*2*2 
supercell of it, and replaced one O atom with my dopant. According to the 
siesta manual, can l use NetCharge and simulatedoping together for my bulk 
structure to simulate doping? Or it is impossible? And if it is incorrect how 
can I dope my bulk structure in siesta?
Yours sincerelyMasoudeh maleki
Dear Masoudeh,
if you explicitly add the dopant atom to your system, you don't need any 
NetChargesince your system is not charged. (The dopant atom comes with its 
correct nuclear chargeand number of electrons). As explained in the manual, the 
NetCharge option is providedto treat charged finite systems (molecules). Best 
regards
Andrei Postnikov



-- 
-- prof. Andrei Postnikov -- tel. +33-372749149 -- cell +33-636565128 --
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
ICPM, 1 Bd Arago - BP 95823, F-57078 Metz Cedex 03, France

Re: [SIESTA-L] Doping

[Andrei Postnikov]

- masoudeh maleki  a écrit :
>Dear sir,
I'm using siesta code to investigate the effect of doping on bandstructure and 
PDOS of periodic crystal of tin oxide with tetragonal unitcell. I made a 2*2*2 
supercell of it, and replaced one O atom with my dopant. According to the 
siesta manual, can l use NetCharge and simulatedoping together for my bulk 
structure to simulate doping? Or it is impossible? And if it is incorrect how 
can I dope my bulk structure in siesta?
Yours sincerelyMasoudeh maleki
Dear Masoudeh,if you explicitly add the dopant atom to your system, you don't 
need any NetChargesince your system is not charged. (The dopant atom comes with 
its correct nuclear chargeand number of electrons). As explained in the manual, 
the NetCharge option is providedto treat charged finite systems (molecules). 
Best regardsAndrei Postnikov

Re: [SIESTA-L] Is it Possible to make a movie in SIESTA?

[Andrei Postnikov]

Dear Rishi,there is my old tool which helps to manipulate the Siesta output 
informationconcerning molecular dynamics or relaxation(available in MD or 
MD_CAR or ANI files) to write more or less directlyinto the animation file 
readable by XCrySDen is available 
fromhttp://www.home.uni-osnabrueck.de/apostnik/download.html( -> Sies2xsf ),or 
in the Siesta distribution, see Util/Contrib/APostnikovSee a short description 
in the README file and at the beginning of the code, md2axsf.fmoreover 
e.g.https://departments.icmab.es/leem/siesta/Documentation/Tutorials/Barcelona-2007/Postnikov-visual.pdffrom
 page 14 on.The XCrySDen nicely allows to run animations and save them into 
animated files;if you don't like XCrySDen or have problems with ityou can 
easily adapt the format of the resulting .axsf file(atomic types and 
coordinates)to any visualising software. Even if you generate a picture "by 
hand" from every selected frame, you'll easily find tools to merge them into a 
movie.Good luckAndrei 
-- prof. Andrei Postnikov -- tel. +33-372749149 ---
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
ICPM, 1 Bd Arago - BP 95823, F-57078 Metz Cedex 03, France



- Rishi Sreedhar <rishisr33d...@gmail.com> a écrit :
>Dear SIESTA community,

I was wondering if it would be possible to get the cartesian output of the 
input molecular cluster whilst SIESTA is running; like say after every 10-20 CG 
optimization moves so that one could gather all those XYZ outputs and treat 
each one as a frame in a movie that visually depicts the system going to it's 
minimum energy configuration starting from the initial state.

‌Does anybody know if such a feature already exists or if it is a tweak that is 
implementable with not much difficulty?
Thank you for your time,
Best Regards,
Rishi Sreedhar.

Re: [SIESTA-L] asterisks and NAN in output

[Andrei Postnikov]

Dear Mogus:I think, the OrderN method is generally much more fragile than 
diagonalisation,so that you should not expect that it goes by itself from 
scratch.In your case, asterisks indicate that your calculation goes astrayearly 
enough (iter 9).In fact you are not obliged to use OrderN just for the sake of 
it,once your system has more than 100 atoms. In many cases, you'llgo ahead just 
fine with systems of many hundreds of atoms anddiagonalisation, and spare 
yourself many troubles. OrderN is justan option to make calculation faster, if 
possible.I'd suggest the following approach:1. Check that the structure is OK. 
It does not have to be pre-relaxed by VASPor whatever, but you should exclude 
error in the units (Bohr / Ang),wrong translation vectors, atoms sitting on top 
of the others, etc.It seems obvious but in fact a large fraction of problems 
comes from his kind of things. Any visualisation of structure would be 
helpful.2. Make a Diagon run and save the density matrix; later on start 
OrderNfrom it, and not from scratch.3. Also from Diagon calculation, check that 
you have a band gap(in the right place), and write down the Fermi energy for 
the use with OrderN.4. Carefully try OrderN to see if it would work fine and 
faster than Diagon.Be prudent with mixing parameters. Look at how the things 
changefrom one iteration to the other and tune the parameters accordingly(to 
have a good yet stable convergence).
Good luckAndrei
-- prof. Andrei Postnikov -- tel. +33-372749149 -- --
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
ICPM, 1 Bd Arago - BP 95823, F-57078 Metz Cedex 03, France


- Mochena, Mogus D. <mogus.moch...@famu.edu> a écrit :

Dear Siesta Users:
>   I am trying to use ORDERN solutionmethod to optimize a quantum dot of CdSe 
> with 90 atoms of Cd and 90 atoms of Se. Actually the system was optimized 
> with VASP and the optimized coordinates were being used to validate SIESTA 
> input parameters prior to
> using SIESTA for larger systems. The diagonal solution method ran fine 
> (single point calculation). When I tried to run the ORDERN solution method, 
> some results (numbers) are getting too large and producing asterisks or a NAN 
> (I think). I don't know how to
> fix the problem. I used the diagonal method results (Md.UseSaveXV) to run the 
> ordern calculation. Please see the output file below. I placed the 
> coordinates at the end of the file.
>
siesta: System type = molecule  

>

> initatomlists: Number of atoms, orbitals, and projectors:180  2520  2880

>

> siesta:  Simulation parameters 
> 

> siesta:

> siesta: The following are some of the parameters of the simulation.

> siesta: A complete list of the parameters used, including default values,

> siesta: can be found in file out.fdf

> siesta:

> redata: Non-Collinear-spin run   = F

> redata: SpinPolarized (Up/Down) run  = F

> redata: Number of spin components= 1

> redata: Long output  = F

> redata: Number of Atomic Species =2

> redata: Charge density info will appear in .RHO file

> redata: Write Mulliken Pop.  = NO

> redata: Mesh Cutoff  =   100.  Ry

> redata: Net charge of the system = 0. |e|

> redata: Max. number of SCF Iter  =  100

> redata: Performing Pulay mixing using= 8 iterations

> redata: Mix DM in first SCF step ?   = F

> redata: Write Pulay info on disk?= F

> redata: Discard 1st Pulay DM after  kick = F

> redata: New DM Mixing Weight = 0.0500

> redata: New DM Occupancy tolerance   = 0.0001

> redata: No kicks to SCF

> redata: DM Mixing Weight for Kicks   = 0.5000

> redata: DM Tolerance for SCF = 0.000100

> redata: Require Energy convergence for SCF = F

> redata: DM Energy tolerance for SCF  = 0.000100 eV

> redata: Require Harris convergence for SCF = F

> redata: DM Harris energy tolerance for SCF = 0.000100 eV

> redata: Using Saved Data (generic)   = F

> redata: Use continuation files for DM= F

> redata: Neglect nonoverlap interactions  = F

> redata: Method of Calculation= Order-N

> redata: Fix the spin of the system   = F

> redata: Maximum number of iterations = 1000

> redata: Relative tolerance   = 0.10D-07

> redata: Eta (Fermi level parameter)  = 0.  Ry

> redata: Radius of LWFs   = 9.5000  Bohr

> redata: Use continuation files for LWF   = F

> redata: Method to build LWFs = kim  

> redata: Dyn

Re: [SIESTA-L] Data in SystemLabel.vectors

[Andrei Postnikov]

- T. Liu <tl...@cam.ac.uk> a écrit :
> Dear All,
> 
> After using Vibra to get the vibrational frequency, I can get such a 
> file called SystemLabel.vectors, and there are many coordinates (or I 
> may call them vectors) in the file as following,
(...)
> Could anyone tell me what these coordinates are about please?


Dear Tao -
Yes, it would make sense to call them (eigen)vectors.
They are relative displacements within the mode in question
times sqrt of mass,
normalised to 1.0
(sum of squares over all atoms and Cartesian coordinates).

Good luck

Andrei Postnikov

> 
> k= 0.000.000.00
> Eigenvector  = 1
> Frequency=-0.108924
> Eigenmode (real part)
>0.6144E-02  0.2053E-01  0.3226E+00
>0.6144E-02  0.2053E-01  0.3226E+00
>0.6269E-02  0.2095E-01  0.3291E+00
>0.6269E-02  0.2095E-01  0.3291E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.4731E-02  0.1581E-01  0.2484E+00
>0.1187E-02  0.3968E-02  0.6235E-01
>0.1187E-02  0.3968E-02  0.6235E-01
>0.1187E-02  0.3968E-02  0.6235E-01
>0.1187E-02  0.3968E-02  0.6235E-01
> Eigenmode (imaginary part)
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00  0.E+00
>0.E+00  0.E+00 -0.E+00
> 
> Could anyone tell me what these coordinates are about please? Are they 
> the deviations (xyz) from the equilibrium position of each atom under 
> this normal mode? What's the Unit for these coordinates?
> Thanks.
> 
> Regards,
> Tao
>

Re: [SIESTA-L] calculate the vibrational frequencies for a molecule with fractional charges

[Andrei Postnikov]

Dear Zhuang,
I'd guess, the issue of charge is technically decoupled from the issue of 
vibrations.
That means, to arrange a calculation for a molecule with fractional
charge is one thing that may pose problems,
but once this is done and the static calculation makes sense
(there is a switch NetCharge;
the documentation says the fractional charges are possible, but warns 
about difficulties),
you go ahead with phonons as you'd do in any other case, with
MD.TypeOfRun = FC. 

Best regards

Andrei

-- prof. Andrei Postnikov -- tel. +33-372749149 -- --
University of Lorraine - Laboratoire de Chimie/Physique - A2MC
ICPM, 1 Bd Arago - BP 95823, F-57078 Metz Cedex 03, France



- Zhuang Bo <zhuangbo2...@outlook.com> a écrit :
> Dear Siesta users,
> 
> I am a  beginer to Siesta and I am trying to use Siesta to calculate the 
> vibrational frequencies for a molecule with fractional charges (for example, 
> 1/2 or 1/3 of charges).
> I found some paper using Siesta to do such kind of calculation, such as 
> 10.1073/pnas.1320210111 and 10.1039/C6RA21476B, but they didn’t provide the 
> procedure in detail that I can follow.
> I have also read through the user manual but couldn’t find the answer.
> 
> Dose anyone here know how to do such kind of calculation using Siesta?
> I really appreciate any help you can provide.
> 
> 
> Zhuang
--

Re: [SIESTA-L] range of the basis function

[Andrei Postnikov]

- toufik esssakhri  a écrit :

> Dear all,> please how to "extract" the range of the basis function in 
> Angstroms from *.psf file?
Dear Toufik,the basis function is NOT in the .psf file (which has to do with 
the pseudopotential);to my knowledge, you can read the basis function extension 
(in Bohr)in the output file, close to the beginning, e.g.:
Vna: Cut-off radius for the neutral-atom potential: 6.351100
Best regardsAndrei Postnikov

Re: [SIESTA-L] Error in Using PAO.SpliNorm

[Andrei Postnikov]

Dear Ajay,
the warning says, certain atoms are too close,
and this most probably happens on applying the translations.
Almost certainly there is an error in structure definition -
I'd suggest to get a view of the atoms in the unit cell,
with and without translations.
SplitNorm is not a problem (at least, not at this stage..) 

Best regards

Andrei Postnikov

- Ajay_Khanna <samthecha...@gmail.com> a écrit :
> Greetings Everyone, I'm wondering if someone can help me in correcting this
> warning
> 
> siesta: WARNING: Atoms 11679 11979 too close: rij =0.402646 Ang
> siesta: WARNING: Atoms 11680 11963 too close: rij =0.516207 Ang
> 
> I'm using PAO.SplitNorm 0.16 and varied it to 0.37
>  But still, I'm getting this warning. Can anybody please tell me what this
> warning about is?? and how can I correct this?
> 
> Thank you.
> Ajay Khanna
> Student

Re: [SIESTA-L] Is relaxation is necessary during k-point optimization

[Andrei Postnikov]

Dear Rajan,as you speak about a ribbon, it looks probable that you'll have 
dispersionin one dimension only, so the k-mesh (to test the convergenceof 
results) will be linear and hardly prohibitively large.I don't understand what 
you mean by "relax my system at each k-point".If you want to test the 
convergence of results as the number of k-pointsincreases, you don't really 
need to relax the structure... Anyway, the result of the convergence test would 
depend on which properties exactlyare you after. To plot energy vs number of 
k-points would hardly be usefulby itself, because it is a function which goes 
down and down and down...You must know where to stop. You'll not be able to 
"evaluate optimized k-point"from the plot.Good luckAndrei Postnikov
- RAJAN SINGH  a écrit :
>Dear Siesta Users> I have a fairly large system of about 50 atoms of graphene 
>ribbon. I want to optimize the same in terms of k-points and mesh cutoff.> But 
>if I relax my system at each k-point, It will cost too much time.> My question 
>is> Can I skip the relaxation fr each k point??> And plot Energies vs kpoint n 
>evaluate optimized k-point??> Please help me in this.With Regarda>Rajan

Re: [SIESTA-L] How to calculate single molecule energies

[Andrei Postnikov]

- I. Camps <ica...@gmail.com> a écrit :
>> Dear Andrei,
>> Thank you very much for your email.
>> What I mean by convergence study is:
> - I choose the accuracy for my energy calculations
> - Do several calculations with different values of meshcutoff energies
> - Plot the system energy (converged within the accuracy selected previously)
> - Select the value of meshcutoff from were the system energy is stable
>> Is this ok?
> []'s,
> Camps

Dear Camps,This is OK up to the last point. What is your definition of 
"stable"?Does the energy have to be stable to within 0.01 eV or 0.1 meV ? This 
would of coursedepend on your problem under study. Moreover, the "system 
energy" as such(the absolute value as it comes out of the calculation) has no 
big meaning;what really interests you is typically the difference of this 
energy from someother energy (that of different phase, different magnetic 
state, etc.).This physically meaningful difference may converge to physically 
meaningful accuracy faster (or, on the contrary, slower)than you might have 
guessed while looking at your [system energy vs cutoff] plot.
Best regardsAndrei 

>
On Mon, Dec 19, 2016 at 8:59 PM, Andrei Postnikov 
<andrei.postni...@univ-lorraine.fr> wrote:
A small comment to the latest remarks,sorry as it must be obvious:it hardly 
makes sense to run the convergence test on a calculation of just one 
molecule(or, in any standalone calculation) because you wont't have a cluewhat 
exactly, and to which accuracy, should be converged.Typically, you are after 
the energy difference between two (or more) cases:atoms vs molecule, free 
molecule + surface vs adsorbed molecule, etc.The  convergence to be checkedis 
that of the meaningful energy difference, to the accuracy you need. 
Best regards
Andrei Postnikov
- I. Camps <ica...@gmail.com> a écrit :
>Dear Riya,

Its always good that, instead selecting those parameters (cutoff energy, cell 
parameters, numbers of k points, etc.) manually, you run convergence test. In 
this way, maybe you get smaller values for them and then less computational 
resources involved.
[]'s,

Camps

>
On Sat, Dec 17, 2016 at 8:23 PM, Riya Rogers <rogers.riy...@gmail.com> wrote:
Dear Dr Camps
>Thank u so much. I m taking very high cutoff greater than 340Ry. N put the 
>Molecule in middle of big lattice.

> My bond lengths are lil away frm exam abt 1.2 to 1.4%.
>Thanking you
>Regards

> Riya
>
On 18-Dec-2016 2:40 am, "I. Camps" <ica...@gmail.com> wrote:


>Can anyone tell me how to calculate energy of one molecule of gas.
>Should I take molecule inside lattice? 
> Or
> Should I follow the following example1:
>http://personales.unican.es/junqueraj/JavierJunquera_files/Metodos/Basic/Basic-exercises.html
SIESTA is designed for periodic calculation, so, you have to put your molecule 
inside a fictional lattice. The lattice vectors should be big enough to not get 
interaction between two neighborhood molecules. Take a look in the second 
example. Also, you can check how SIESTA is treating your system (molecular, 
surface, slab, etc.) in the output file.

>Does k points n meshcutoff and basis set info req?Yes.

K vectors has meaning only for periodic system, so, in case of a molecule, only 
one is needed (the gamma point (0,0,0).
For the meshcutoff, you should do the convergence study as usually.
The basis set should be selected according the accuracy you want or need (as in 
any SIESTA calculation).
Regards,
-- 
[]`s
>Camps
>
>
>

Re: [SIESTA-L] How to calculate single molecule energies

[Andrei Postnikov]

A small comment to the latest remarks,sorry as it must be obvious:it hardly 
makes sense to run the convergence test on a calculation of just one 
molecule(or, in any standalone calculation) because you wont't have a cluewhat 
exactly, and to which accuracy, should be converged.Typically, you are after 
the energy difference between two (or more) cases: atoms vs molecule, free 
molecule + surface vs adsorbed molecule, etc. The  convergence to be checkedis 
that of the meaningful energy difference, to the accuracy you need. Best 
regardsAndrei Postnikov
- I. Camps  a écrit :
>Dear Riya,

Its always good that, instead selecting those parameters (cutoff energy, cell 
parameters, numbers of k points, etc.) manually, you run convergence test. In 
this way, maybe you get smaller values for them and then less computational 
resources involved.
[]'s,

Camps
 
>
On Sat, Dec 17, 2016 at 8:23 PM, Riya Rogers  wrote:
Dear Dr Camps
>Thank u so much. I m taking very high cutoff greater than 340Ry. N put the 
>Molecule in middle of big lattice.

> My bond lengths are lil away frm exam abt 1.2 to 1.4%.
>Thanking you
>Regards

> Riya
>
On 18-Dec-2016 2:40 am, "I. Camps"  wrote:


>Can anyone tell me how to calculate energy of one molecule of gas.
>Should I take molecule inside lattice? 
> Or
> Should I follow the following example1:
>http://personales.unican.es/junqueraj/JavierJunquera_files/Metodos/Basic/Basic-exercises.html
SIESTA is designed for periodic calculation, so, you have to put your molecule 
inside a fictional lattice. The lattice vectors should be big enough to not get 
interaction between two neighborhood molecules. Take a look in the second 
example. Also, you can check how SIESTA is treating your system (molecular, 
surface, slab, etc.) in the output file.

>Does k points n meshcutoff and basis set info req?Yes.

K vectors has meaning only for periodic system, so, in case of a molecule, only 
one is needed (the gamma point (0,0,0).
For the meshcutoff, you should do the convergence study as usually.
The basis set should be selected according the accuracy you want or need (as in 
any SIESTA calculation).
Regards,
-- 
[]`s
>Camps
>
>
>

Re: [SIESTA-L] siesta

[Andrei Postnikov]

Dear 木材,you might find useful to have a look 
athttp://departments.icmab.es/leem/siesta/Documentation/Tutorials/Lyon-2007/Postnikov-visual.pdf(from
 page 9 on),best regardsAndrei Postnikov
- 木材  a écrit :
>Dear all When complete the calculation of siesta and get a file 
>selected.WFSX,I know it is for LUMO and HOMO,but how to deal with it and plot? 
>I really hope to get your hlep,hope your comment.

Re: [SIESTA-L] WRONG LOCAL MAGNETIC MOMENT OF Fe IN FeO.

[Andrei Postnikov]

Dear Anju,1) you lattice parameter seems to be too small; you define 3.04 
Angwhereas in wüstite it seems to be around 4.3 Ang. Hopefully this is the 
source of your problem.2) Before playing around with U (up to U=14 eV, 
indeed!), why don't you check your calculation against the known LDA results? 
(Not only the magnetic moment but the placement of bands...)3) When asking for 
advise in searching for an error, it might be a good idea to provide an output 
along with your input file... Good luckAndrei Postnikov
- Anju Saroha  a écrit :
>Hi all,

I have tried siesta-trunk-587 (latest version) but I am getting the same error. 
In the FeO calculation, simply to check, I gave ferromagnetic arrangement of 
spins for Fe. The expected local magnetic moments of Fe should come out to be 4 
Bohr magneton but I am getting 0.001 Bohr magneton. Even for U=6 eV there is no 
band-gap. 

Can anyone please help me in resolving this issue? I have tried my best but 
didn't understand where is it wrong. I am looking forward for some suggestions. 
Please help me out. Here is the attached .fdf file.

On Fri, Nov 25, 2016 at 9:28 PM, Anju Saroha  wrote:
Sorry I forgot to attach the input file. I am attaching the FeO.fdf file 
herewith. Please find the attachment.



On Fri, Nov 25, 2016 at 9:11 PM, Anju Saroha  wrote:
Hi all,

I am using siesta-trunk-527. I did a test run on FeO to check the installation. 
But I have encountered an error. Though FeO is an antiferromagnetic material, I 
gave the calculation for ferromagnetic case. The net magnetic moment of the 
unit cell is expected to be 8 Bohr magneton but I have got it as 0.5 Bohr 
magneton.

The local magnetic moments of Fe are coming around 0.001 Bohr magneton 
(expected are 4 Bohr magneton). 

Has anyone faced a similar problem? Is it a bug in the code? If not, then what 
can lead to such an erroneous result? Even FeO is coming out to be metallic 
till U = 14 eV. Though it should become an insulator for U to be around 6 eV. 
Can anyone please help me to understand the error? I tried a lot to solve it 
but I am now clueless. Any help in the direction will be appreciated. I am 
curiously looking forward for suggestions. Thanks in advance.

Regards

-- 
Anju Saroha
Research scholar
IIT Madras
Chennai-600036
India.

Re: [SIESTA-L] How to extract the force constant after siesta'scalculation for thermal conductance calculation by using nonequilibriumGreen's method

[Andrei Postnikov]

Dear Benhu Zhou,
you want a lecture on how phonon calculation in Siesta works,that will replace 
all previous documents and tutorials. OK, I'll try to make it short:1. You have 
unit cell that contains M atoms. For graphene, silicon, NaCl M=2.2. If you want 
phonon dispersions, you construct a supercell that is a replicatedunit cell (in 
one dimension, or two, or three, depending on which dispersioninterests you and 
how far-fetching the force constants are). The supercell containssay N atoms 
(an integer multiple of M). If you want only Gamma phonon and no dispersion, 
supercell = unit cell and M=N.3. The FC calculation is organized as several 
inserted loops.The outer loop is over atoms in the unit cell, 1 to M. (You can 
performin your calculation any part of this loop, from atom FCfirst to FClast,1 
<= FCfirst <= FClast <= M, and then concatenate them by hand). Every run of 
this loop adds 6*N lines to the .FC file.The next inner loop is over 
displacements, 1 to 6.For every displacement, the forces are calculated on all 
atoms in the supercelland written to the .FC file as a block of N lines, 3 
columns (for Cartesian componentsof forces) each.As the FC calculation runs, it 
explains in the output file what it is doing right now,and you can inspect how 
the blocks are added to the.FC file.Good luckAndrei
- 周本良 <blz...@hunnu.edu.cn> a écrit :
>Dear Andrei Postnikov,
>   Thank you very much for your answer and offering relevant materials to me!
>Yes, your offer materials has bee read in detail, and the similar issue 
> was also introdeced by Bedamani (see the website 
> https://www.mail-archive.com/siesta-l@uam.es/msg04969.html), and a reply was 
> given by Roberto.  Part explanations are follows:
> "Let's say FC.fdf file built by fcbuild contains N atoms, then each of them 
> gets a block of 6N rows by 3 columns in the *.FC file, representing the 
> forces on every atom, whenever every  atom of your basic, reference, unit 
> cell moves a small distance  (delta) along +x,-x,+y,-y,+z,-z. The numbers 
> written are not the forces themselves but forces/delta. Thus, if your unit 
> cell contains M atoms, the *.FC file will
 comprise M*6N lines."
>Blame on me for being too stupid, I still have some question.
>(1) as mentioned above, "a atom have 6 rows by 3 columns in the *.FC file, 
>representing the forces on this atom, whenever the atom of your basic, 
>reference, unit cell moves a small distance  (delta) along +x,-x,+y,-y,+z,-z." 
>Dose the 6 rows represent the force constant for a atom in the unit cell when 
>the atom moves a small distance along +x,-x,+y,-y,+z,-z?  if it is true, what 
>the mean for the first column, the second column and the third column are? 
>what the mean for the words "basis" and "reference" are?
>(2) "if unit cell contains M atoms, the *.FC file will comprise M*6N lines." I 
>don't undstand why the total lines are M*6N, I think the total lines are 6*N 
>because N is the total number of atoms in the supercell.
>
> 
>Please help me, thanks!
>Best regards,
>Benhu Zhou
>
>-- Original --
>
>From:  "Andrei Postnikov"<andrei.postni...@univ-lorraine.fr>;
>Date:  Sat, Aug 20, 2016 09:59 AM
>To:  "siesta-l"<siesta-l@uam.es>;
>

>Subject:  Re: [SIESTA-L] How to extract the force constant after 
>siesta'scalculation for thermal conductance calculation by using 
>nonequilibriumGreen's method
> 
>
> 
>Dear Benhu Zhou,

>for the answer to (1), you can consult page 15 of

>www.home.uni-osnabrueck.de/apostnik/Lectures/APostnikov-Phonons.pdf
>(also, I think the issue was discussed many enough times in the mailing list).

>

>Best regards

>

>Andrei Postnikov

>

>- 周本良 <blz...@hunnu.edu.cn> a écrit :
>
>Dear all,
>
>
> I am trying to use the force constant output from siesta for a phonon 
> transmission calculation for a graphene nanoribbon by using non-equilibrium 
> Green's function. The system for a graphene nanoribbons is consists of a 
> graphene nanoribbon (a unit cell have 8 atoms ) connected to two graphene 
> leads. Now I obtain the file graphene.Fc (24 atoms, 144 rows, 3 column, and 
> periodic boundary conditon is applied in the heat transport direction) with 
> the force constat matrix. But I have some questions about the file 
> graphene.FC.
>
>
>(1)what is the mean for the elements in the graphene.FC file? For example, 
>what is the mean for the first column or the sencond culumn or the third 
>column?  and what is the mean for the first row or the second row and so on?  
>I have read the guide in detail but I still do not understand. 
>
>
>(2) How I can transform the force constant to a  matrix with 24X24 (basicaly

Re: [SIESTA-L] H2S Ground state has imaginary frequencies

[Andrei Postnikov]

Dear Anja,
the problem of negative frequencies cannot be helped by argumentationabour 
agreement with experiment or other calculations. In fact it is a technical 
problem which you can explore very precisely,whatever your pseudos or bases. 
You have a certain total energy surface(as function of all coordinates). You 
apply a CG algorithm and arriveat what seems to be a minimum. The negative 
frequency hintsthat the starting point was in fact NOT a minimum, and 
indicateswhere to search better. This strategy is fool-proof, unless:i) your 
total energy surface is pathological, for whatever reason;ii) the numerical 
noise (fluctuations of total energy) in the course of relaxationis comparable 
with the energy differences you are after.The case (ii) can be systematically 
cured by increasing the MeshCutoff.Now as I look at your translation vectors 
(1.803645   15.6477590.331091), (11.8583939.3745580.435976),  
(0.5925850.553982   24.937136)it sems that you let them vary in the course 
of relaxation? -(Hopefully NOT in the course of phonon calculations? - Sorry).I 
see that your simulation boxis quite distorted - not an error by itself but I'd 
quess that your initialchoice was more "rectangular". A considerable 
deformation of the boxfor me is an indication that theMeshCutoff is too 
small.The fact that you need to go to unreasonably high MeshCutoffmight be a 
consequence of unfortunate choice of pseudos etc.but in principle you can 
eliminate the negative frequency problemjust by increasing the mesh, whatever 
the other settings.Best regardsAndrei - Anja Foerster <foers...@msu.edu> a 
écrit :

> Dear Andrei Postnikov,
> Thank you for your reply. I already followed the movement of the imaginary 
> frequency and used the new structure as the input for a new 
> optimization/frequency run. Essentially, I still end up with the same 
> optimized structure and frequency. If I compare my structure with the 
> literature values for the angle and distance, it is in good agreement for a 
> DFT calculation and as the by hand change of the angle and distance show, it 
> is the energetic minimum.
> As for the suggestion to increase the size of the box: While for the simple 
> H2S molecule I can increase the box size to test if dipole-dipole 
> interactions are the cause, the H2S molecule is just a test system. What I 
> really want to calculate is H2S reacting with a 2-dimensional surface (about 
> 50 atoms). Here, I am not able to increase the size of the box due to the 
> considerable increase of the atom number and the calculation time in such a 
> case, especially for the frequency run. Therefore, I would be grateful if 
> there was a way to get the right frequencies for H2S while keeping the box 
> size of the 2D surface, which should be sufficiently large for H2S
> the box size of the 2D surface in Ang : > 1.803645   15.6477590.331091
> 11.8583939.3745580.435976
> 0.5925850.553982   24.937136
> Thank you,

> Anja

> Quoting Andrei Postnikov <andrei.postni...@univ-lorraine.fr>:
> > Dear Anja,
> > in my experience, large negative frequencies are either trivial
> > "noise"(too low mesh cutoff - apparently not your case), or, more
> > frequently, a hint that you are not really at equilibrium. A good
> > idea would be to look at the eigenvectorof the "most negative" mode,
> > extract the relative displacements(not to forget the sqrt of masses)
> > and give your atoms a small kickin this direction before resuming the
> > relaxation. The negative (imaginary)frequency means - you displace
> > the atoms in this mode, and the energy lowers.This might help you to
> > find a better energy minimum,calculate there the frequencies again,
> > etc.Another trouble (or, the reason for the behavior described
> > above)could be the dipole-dipole interactions between your
> > moleculesif the box is not large enough. However, the d-d interaction
> > is long-ranged, so you'll have some amount of trouble at any box
> > size.I don't remember if there are efficient switches to cancel it in
> > this context.For a molecule like H2S, you would normally get 3 rigid
> > displacementand 3 rigid rotation modes, all at zero frequency. It is
> > relatively easyto bring the displacement modes to zero (by
> > brute-force increasingthe mesh cutoff), whereas bringing the rotation
> > modes to zero can bea hell, apparently because the "Siesta space" is
> > not isotrope enough (interaction with replicas of the molecule, the
> > rigid structure of mesh, ..?)Best regardsAndrei Postnikov

> > - Anja Foerster <foers...@msu.edu> a écrit :

> >

> >> Dear all,

> >> I performed a geometry optimization of H2S and then used the Vib

Re: [SIESTA-L] H2S Ground state has imaginary frequencies

[Andrei Postnikov]

Dear Anja,
in my experience, large negative frequencies are either trivial "noise"(too low 
mesh cutoff - apparently not your case), or, more frequently, a hint that you 
are not really at equilibrium. A good idea would be to look at the 
eigenvectorof the "most negative" mode, extract the relative displacements(not 
to forget the sqrt of masses) and give your atoms a small kickin this direction 
before resuming the relaxation. The negative (imaginary)frequency means - you 
displace the atoms in this mode, and the energy lowers.This might help you to 
find a better energy minimum,calculate there the frequencies again, etc.Another 
trouble (or, the reason for the behavior described above)could be the 
dipole-dipole interactions between your moleculesif the box is not large 
enough. However, the d-d interaction is long-ranged, so you'll have some amount 
of trouble at any box size.I don't remember if there are efficient switches to 
cancel it in this context.For a molecule like H2S, you would normally get 3 
rigid displacementand 3 rigid rotation modes, all at zero frequency. It is 
relatively easyto bring the displacement modes to zero (by brute-force 
increasingthe mesh cutoff), whereas bringing the rotation modes to zero can bea 
hell, apparently because the "Siesta space" is not isotrope enough (interaction 
with replicas of the molecule, the rigid structure of mesh, ..?)Best 
regardsAndrei Postnikov
- Anja Foerster  a écrit :

>Dear all,
> I performed a geometry optimization of H2S and then used the Vibra tool and a 
> FC run to calculate the frequency. For the optimized structure I receive 5 
> negative (imaginary) frequencies (-228.9332   -135.1068-46.3373 
> -0.1270 -0.0328). Since it is a small molecule, I varied both the 
> H-S-distance and the H-S-H angle. As you can see in the attached pdf, 1.38 
> Ang is the H-S distance that corresponds to the energetic minimum. However, 
> only when the H-S-distance is increased to 1.40Ang are the frequencies 
> getting close enough to 0 for the structure to be labeled as "ground state" 
> based on the frequency results.
> I already tried to increase the meshcutoff value to 800, which did not solve 
> the problem. My other critical values: MD.MaxForceTol is  0.001 eV/Ang and 
> DM.Energy.Tolerance  is 1.0d-5 eV. I set the MD.MaxForceTol to a stricter 
> 0.0001 eV/Ang in the geometry optimization calculation, however even after 
> 400 steps is not converging. For MD.MaxForceTol 0.001 eV/Ang I receive “Tot   
>  0.0026530.001420   -0.001683” in the siesta output for “ siesta: Atomic 
> forces (eV/Ang)”. So the forces should be okay for a frequency run.
> Does anyone have an explanation why the energetically lowest state has 5 
> imaginary frequencies whereas an energetically higher geometry has 
> effectively no imaginary frequencies?
> Thank you for your help,
> Anja
> P.S.: I'm using siesta version 3.1

>

Re: [SIESTA-L] Cholesky factorisation in rdiag error due to periodicity

[Andrei Postnikov]

Hi Arturo,
you definedAtomicCoordinatesFormat Fractionalso that your atoms Nr 1 and 2
%block AtomicCoordinatesAndAtomicSpecies
4.00 0.00 1.00 3
3.00 0.00 1.00 3
as well as many others, are at identical positions. You know, in "fractional" 
arithmetics 4=3 (=2=1=0)You can see this explicitly as the coordinates are 
printed in Bohr: 
siesta: Atomic coordinates (Bohr) and species
siesta: 152.22128 0.0 10.12043 3 1
siesta: 114.16596 0.0 10.12043 3 2
- these coordinates are separated by exactly the 1st translation vector(20.138 
Ang).Indeed, this is because of periodicity (which is implicit in 
Siesta),however this is not an error but exactly the structure as you decribed 
it...The first step in solving the problem would be to figure outwhich geometry 
you really want to use (how long the translation vectorsand where are the atoms 
in the cell: placed from one another at distance =1.0 in which units?) Good 
luckAndrei Postnikov- Arturo Toro  a écrit :
>Hi all,

While attempting to run a molecular dynamics simulation of an interface between 
a hexagonal system and a cubic system, I received the above mentioned error 
message. I then looked through the .out file and saw that there were numerous 
atoms which, although visibly located at different atoms, are said to be 0 
angstroms apart.

I believe this to be the result of an innate periodicity function within SIESTA.
Is there a way to avoid this error besides manually attempting to delete all 
the conflicting atoms?

Thanks again,
Arturo Toro

Re: [SIESTA-L] bulk modulus by Postnikov's tool: Murnaghan equation of state fit.

[Andrei Postnikov]

Dear Manjeet:
the "concept" is to add one line of code that calculatesvolume (in Bohr^3) from 
your lattice constant (in Ang).There is an example (commented) for bcc 
lattice.When programming this transformation, take care that the volumeyou 
calculate corresponds to the energy you have.Good luckAndrei Postnikov
- Manjeet Bhatia <manjeetbhati...@gmail.com> a écrit :
>Dear Andrei Sir,

Thank you for your explanation. Yes, my input is lattice constant in ang and 
energy in eV. This means that there is no concept of lattice constant as input 
in murnagan fit, only volume and energy in appropriate units will give bulk 
modulus.


Thank you.  

On Thu, Aug 11, 2016 at 9:14 PM, Andrei Postnikov 
<andrei.postni...@univ-lorraine.fr> wrote:
Dear Manjeet:
The mur_fit script expects the input data to be passed
in two columns:
(1) volume in Bohr^3, (2) Energy in Ry.
These are the units internally used, in which the results
are provided. In you case: 
#  Equilibrium volume, a.u.^3 :   31.125644
#  Minimum energy, Ry : -614.470876Now, it might be in the user's 
convenience to pass instead
the volume and/or energy in other units, and convert them
into Bohr^3 / Ry once they are read in. To this end, some"templates" are 
provided in the source file,(before the i=i+1 statement)
i.e. conversion operations
which can be (un)commented. In the mur_fit.f you attached,
the volume is presumed to be in Ang^3 and the energy in eV.
Indeed, your data file has the volume corresponding to
the lowest energy = 4.61 (presumably Ang^3),
and the corresponding energy = -8362.950245 (presumably eV).
You'll find nearly these values, expressed in Bohr^3 and Ry,in the output from 
your test. The bulk modulus is crazy (and does not fitinto the print format) 
because the units (at least in the 1st column)are most probably wrong: for me, 
the values like 4.61 do not look like the volume in Ang^3.It may well be the 
lattice constant, in which case you'd need to activate a line to convert it 
into volume in Bohr^3. 

Summarizing:
1. Use the units of your convenience, but check that the results
after the fit are correct.
2. Make sure that the energy (as in column 2) is given 
exactly per volume (as in column 1).Do not forget that e.g. the fcc cube has 4 
atoms
whereas its primitive cell has only one etc.

Good luck

Andrei Postnikov



- Manjeet Bhatia <manjeetbhati...@gmail.com> a écrit :
>Thank you for your reply,please find attached files.

On Thu, Aug 11, 2016 at 6:20 PM, Anita Rani <ranianit...@gmail.com> wrote:
Please send me your murn.f file and .DAT file.

> I will check and then suggest the changes.
>
On Aug 11, 2016 12:20, "Manjeet Bhatia" <manjeetbhati...@gmail.com> wrote:
Dear Siesta users,

I am getting volume, energy and pressure in corresponding units but not getting 
bulk modulus, please suggest if i have to change something in mur_fit.f fortran 
file. 
 
My values are:


#  Equilibrium volume, a.u.^3 :   31.125644
#  Minimum energy, Ry : -614.470876
#  Bulk modulus, Kbar :
#  BP :7.709490




>
>

Re: [SIESTA-L] deal with the .RHO file of siesta

[Andrei Postnikov]

Dear zyynihao,
rho2xsf is not built from a single source file rho2xsf.fbut demands some 
others, that's why a Makefile is providedin the directory which contains the 
source. Chances are,the Makefile will catch up the compiler settings from the 
Siestainstallation, or otherwise you define FC by hand, e.g.,FC = ifortFFLAGS = 
-Oand type "make". It is not necessary to install and run the scriptson the 
same platform where you run Siesta; you can transfer the .XV, .RHO etc. filesto 
your local computer and run rho2xv there (if the coding of binary files is 
compatible,of course).Another hint: there is a README file in the 
distribution.In case of doubt / problems, go get the code from 
http://www.home.uni-osnabrueck.de/apostnik/download.htmlGood luckAndrei 
Postnikovan you does not compuile
- 木材  a écrit :
>Dear all !
>when open Utile/contribe/Apostnikov and find the rho2xsf.f ,then I carry out 
>the order "ifort -o rho2xsf  rho2xsf.f" and it show in the picture,I hope for 
>geting your hlep,how should I do?

Re: [SIESTA-L] Question: Proper SolutionMethod for relaxing large systems

[Andrei Postnikov]

Dear Bahareh:> But I recieve error at the very beginnig of the calculation -  
what kind of error? (What was the error message)?> and the job terminates 
badly!- how badly?It is difficult to help you as you give no details, but let 
me guess. You are using some old version of Siestawhere the SolutionMethod by 
default was diagon for less than 100 atomsand OrderN for beyond 100 atoms. 
That's OK, just define SolutionMethod diagonand go ahead. Yes you can. If the 
calculation is too heavy that's another story;you can always try OrderN later 
(if your system has a gap, etc. etc.)Good luckAndrei Postnikov

- Bahareh Fakhraei <fakhraei.baha...@gmail.com> a écrit :
> Thak you so much dear Andria,
So you suggest that I can use diagon for a system consisted of about 730 
atoms?But I recieve error at the very beginnig of the calculation and the job 
terminates badly!I thought it might be because of that. But if that's ok with 
diagon, I must go over it again and give it another shot.
Sincerely yours,
B. Fakhraei

On Tuesday, August 2, 2016, Andrei Postnikov 
<andrei.postni...@univ-lorraine.fr> wrote:
Hi,

> I think there must be some misunderstanding:

> the Siesta rule of thumb about ~100 atoms implies only that

> for large enough system, O(N) might ultimately provide

> better scaling than diagon, and ~100 is a rough estimate

> of when the crossing of ~N^3 and ~N tendencies happens in practice.

> However, O(N) only works with large enough band gap, and is otherwise fragile.

> There is nothing wrong about using diagon to systems of any size,

> and 100 atoms are not that many to make it markedly unpractical.

> I'd say, go ahead with diagon as long as the calculation time

> is not a problem. If it starts to become an issue

> and you think you can gain by switching to O(N),

> do some test on whether it works correctly and is indeed faster.

>

> Best regards

>

> Andrei

>

>

> - RCP <pasia...@cnea.gov.ar> a écrit :

> > Hi,

> >

> > I'm not an expert here, but the naive answer is "use order N

> > method". However, that depends on the physics of your system,

> > for instance, O(N) cannot be applied to metals.

> >

> > Regards,

> >

> > Roberto

> >

> >

> > On 08/01/2016 05:42 PM, Bahareh Fakhraei wrote:

> > > Hello every one,

> > >

> > > For relaxing the large systems (more than 100 atoms) in SIESTA, seemingly

> > > *diagon* is not a good option for solution method. Why and What should be

> > > used instead?

> > >

> > > yours,

> > >

> > > B. Fakhraei

> > >

> > > A beginner to SIESTA program

>

> --

>

Re: [SIESTA-L] Question: Proper SolutionMethod for relaxing large systems

[Andrei Postnikov]

Hi,
I think there must be some misunderstanding:
the Siesta rule of thumb about ~100 atoms implies only that
for large enough system, O(N) might ultimately provide
better scaling than diagon, and ~100 is a rough estimate
of when the crossing of ~N^3 and ~N tendencies happens in practice.
However, O(N) only works with large enough band gap, and is otherwise fragile.
There is nothing wrong about using diagon to systems of any size,
and 100 atoms are not that many to make it markedly unpractical.
I'd say, go ahead with diagon as long as the calculation time
is not a problem. If it starts to become an issue 
and you think you can gain by switching to O(N),
do some test on whether it works correctly and is indeed faster.

Best regards

Andrei
 

- RCP  a écrit :
> Hi,
> 
> I'm not an expert here, but the naive answer is "use order N
> method". However, that depends on the physics of your system,
> for instance, O(N) cannot be applied to metals.
> 
> Regards,
> 
> Roberto
> 
> 
> On 08/01/2016 05:42 PM, Bahareh Fakhraei wrote:
> > Hello every one,
> >
> > For relaxing the large systems (more than 100 atoms) in SIESTA, seemingly
> > *diagon* is not a good option for solution method. Why and What should be
> > used instead?
> >
> > yours,
> >
> > B. Fakhraei
> >
> > A beginner to SIESTA program

-- 

Re: [SIESTA-L] graphene

[Andrei Postnikov]

In graphene, the bands cross the Fermi energy in the K point. If your k-point 
sampling is sparsebut it so happens that it includes exactly the K point,you'll 
have a peak at the Fermi energy (as in your case).As either the density of 
k-point would increase, orthe tetrahedron integration done instead of 
sampling(not implemented in Siesta, to my knowledge),the importance of the K 
point will diminish,and the "correct" structure of DOS (smoothly converging to 
zeroat the zero gap) would gradually manifest itself.
- Zahra Talebi <talebi_z2...@yahoo.com> a écrit :
>thanks for your answer, 
but you know graphene is zero band gap material, for graphene with 200 atom DOS 
diagram shows the zero band gap but for 288 atom my DOs diagram show a pick on 
the point of fermi energy. 
do you think that it will correct if I increase the kpoints. 
this is my question



 On Monday, July 25, 2016 8:29 PM, Andrei Postnikov 
<andrei.postni...@univ-lorraine.fr> wrote:


Hi Zahrayour both DOS distributions are essentially quite close,only that one 
of them (that with 288 atoms) uses much less k-points(the 6*6 mesh). The 
resulting energy values are too sparse,that's why you see strongly fluctuating 
DOS curve.The other calculation (with 200 atoms) uses much more k-points (10*10 
mesh)and the resulting DOS is much more smooth.Solution:- increase the density 
of k-points (at least just for PDOS calculation, see the switchPDOS.kgrid 
Monkhorst Pack )and/or- increase the broadening of DOS by e.g.setting 
Electronic Temperature, or set the appropriate broadening parameter 
inProjectedDensityOfStatesblock.The exact "fix" depends on what exactly do you 
want:an accurate smooth curve or merely a reasonably looking one,without strong 
unphysical fluctuations.After all, this is just about the graphic 
representation...Best regardsAndrei- Zahra Talebi <talebi_z2...@yahoo.com> 
a écrit :>hi every bodyi am working on graphene surface, and I got a problem. I 
am trying the different super cell. for supercell with 200 atom I got the right 
DOS and band file. but for supercell with 288 atom I got wrong DOS but right 
band. I am attaching the fdf file to my e-mail. does any body knows the 
problem. I also attached the DOS file I really appreciate your guidezahra

Re: [SIESTA-L] graphene

[Andrei Postnikov]

Hi Zahra

your both DOS distributions are essentially quite close,
only that one of them (that with 288 atoms) uses much less k-points
(the 6*6 mesh). The resulting energy values are too sparse,
that's why you see strongly fluctuating DOS curve.
The other calculation (with 200 atoms) uses much more k-points (10*10 mesh)
and the resulting DOS is much more smooth.

Solution:
- increase the density of k-points (at least 
just for PDOS calculation, see the switch
PDOS.kgrid Monkhorst Pack )
and/or
- increase the broadening of DOS by e.g.
setting Electronic Temperature, or 
set the appropriate broadening parameter in
ProjectedDensityOfStates
block.

The exact "fix" depends on what exactly do you want:
an accurate smooth curve or merely a reasonably looking one,
without strong unphysical fluctuations.
After all, this is just about the graphic representation...

Best regardsAndrei
- Zahra Talebi  a écrit :
>hi every body
i am working on graphene surface, and I got a problem. I am trying the 
different super cell. for supercell with 200 atom I got the right DOS and band 
file. but for supercell with 288 atom I got wrong DOS but right band. I am 
attaching the fdf file to my e-mail. does any body knows the problem. I also 
attached the DOS file I really appreciate your guide
zahra

[SIESTA-L] fdf parameters

[Andrei Postnikov]

Dear Amrish:
Let's sort things out:1. The energy shift is about constructing your basis 
functions. There is nothing you need to teston a huge supercell. Make the usual 
tests on small (minimal; single cell) systemsto be sure that your basis set 
construction permits you to obtain the reasonable description(of whatever your 
aims are: lattice parameters, bulk modulus, ...)2. The mesh cutoff, as it seems 
to me, a priori is not supposed to depend on the supercell size.Again, after 
you test it on prototype small systems, you can use the same valuefor large 
supercell.3 my workstation is not so powerful to deal with system like 500 
atoms... - then:why use such workstation? Or: why use 500 atoms? Or: why use 
Siesta?Best regardsAndrei

i am facing problem in my fdf file to take value of some parameters like mesh 
cutt off and energy shift parameters. my workstation is not so powerful to deal 
with system like 500 atoms so i don't want to waste time while test 
calculations on these parameters.can any one explain an alternate way to get 
approximate but the accurate value of mesh cutt off and energy shift for system 
consisting 500 atoms.
-- 
warm regards
amrish sharmaresearch scholarDr. Isha Mudahar grouppunjabi university 
patialahii all 

Re: [SIESTA-L] script for energy vs lattice constant convergence.

[Andrei Postnikov]

Dear Javier, Mayuri:
There was a comment on the related issue
(how to run a sequence of similar calculations, NOT how to plot a curve)
some time back, please see
http://www.mail-archive.com/siesta-l%40listserv.uam.es/msg02458.html
with the batch script integrated into the text of that message 
(from 
#!/bin/sh
to
echo $0" "$1" finished on "$dte >> $wrk0/LOG
)

The script is still operational...
Best regards
Andrei


On Wednesday, July 20, 2016 19:04 CEST, Javier Junquera 
 wrote:
 
Dear Mayuri:

If you want to run various repetitive calculations changing the lattice 
constant in a compact form,
you can implement it yourself starting from some of the examples provided in 
the Util/Scripting directory in the Siesta package.

You can also take a look to the instructions and examples provided in 
Util/JobList.
It contains utilities to organize, to run or queue,
and to collect results of multiple siesta jobs. It is particularly suited
for preliminary convergence tests, and it uses only fortran codes and
minimal shell scripts.

Once you have the energy versus lattice constant curve, you can use the script 
to fit it to the Murnaghan equation of state.
You can download it from:

http://personales.unican.es/junqueraj/
and follow the links:

A self-explained Siesta tutorial
Set of self-explained Siesta exercises
Computing structural and electronic properties of materials
Lattice constant, bulk modulus and equilibrium energy of solids.

Hope this helps,

Javier





> hi all siesta users,
> i want energy vs lattice constant graph. so i need script which is compatible 
> with siesta, i search about how to make script to get minima in ene vs Latt. 
> Cons. but proper result is not obtained.
> please help me..