Re: [SIESTA-L] calculation w/h transition metal
On Sat, 11 Oct 2008, HeeSung Choi wrote: | Dear all, | | Is it possible to calculate a structure wiht transition metal in SIESTA? Dear HeeSung Choi: Yes it is | As I know, transition metal calculation in SIESTA is impossible, You know wrong; there is an impressive number of counter-examples. Check Publications on the Siesta homepage. | because it can not be converged. Yes it can; but may need care. The care regards the fact that in order to converge (any) metal you'd need a quite dense k-mesh, very different from a default value. | My structure is a simple tetrahedron structure with tungsten. OK | Also how do I make Pseudo-potential file for transition metal? Is it | different from other material? Check Pseudos Bases section of Siesta homepage, to begin with. | If anyone did like a this calculation, please give me a help. Which this? Check e.g. http://arxiv.org/abs/cond-mat/0109540 for the calculation of metallic Fe with Siesta, convergence issues etc. For tungsten (in other relation, for a metalloorganic systems) I ended up with a pseudopot which includes 4f in the valence states, which might be a not so obvious choice, but otherwise the pseudopot consruction is pretty straightforward. Good luck Sincerely Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Chimie, Physique et Mat\'eriaux, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] How to define Bulk Modulus
Nidhi: if you are only interested in the E(a) curve, why are you doing MD annealing? You only need a sequence of ~10 static calculations MD.TypeOfRun CG# Default value MD.NumCGsteps 0 (the TypeOfRun does not really matter, but you do not need more than the initial step in whatever kind of MD!) for values of lattice constants within several % around the equilibrium, then you fit the results (total energy) to the Murnaghan (or other) equation of state and you get bulk modulus and its derivative immediately from the fit. If you don't have your favourite fitting progran use mine, http://www.home.uni-osnabrueck.de/apostnik/Software/mur_fit.f (have a look at the fitting curve vs. your data points before believing in the numbers you get from the fit!) The interval of A= 5 to 7 Ang might be OK to localize the minimum, but for a clean fit you'd need much more densely situated points just around the minimum. Good luck, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Chimie, Physique et Mat\'eriaux, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+ On Wed, 23 Apr 2008, Nidhi Sharma wrote: | Dear users, | | I am trying to find the Energy vs alat curve for the B1 phase of LaAs. For this I set the MD. type of run Anneal and Target pressure set to 0.0 GPa and vary the value of alat from 5 to 7 Ang in steps of 0.1. | | MD.TypeOfRun Anneal | MD.AnnealOption Pressure | MD.TargetPressure 0.0 GPa | MD.TauRelax 100.0 fs | MD.BulkModulus100.0 Ry/Bohr**3 | MD.MaxStressTol 0.1 GPa | MD.VariableCell .true. | | I have given default value of Bulk Modulus in .fdf file. What values should I given for MD.BulkModulus in .fdf file. | I got the different energies and pressure for different alat and got the optimal alat. The .out file also report some data in Ry/Bohr**3. I think it is bulk modulus and how we will get the first pressure derivative of Bulk Modulus ? | | Please guide me. | | Thanks in anticipation, | | Nidhi Sharma |
Re: [SIESTA-L] Question about spin configuration
Dear Cheng-Pin Yen (just for keeping our mailing list browsable: you answer to the list, as we were discussing the problem of setting spin at different orbitals in a desired manner, and using DMtune to this end:-) On Tue, 1 Apr 2008, Cheng-Pin Yen wrote: | Thanks a lot for your helping. | I will tune InitSpin block first in fdf input file. If it solves your problem, that's the best solution. | DMtune tools seems to be difficult to use for me ... | I have compiler DMtune in linux . | But I do not know what option to choose for my system . | Is option 5 or 7 right way to do this ? In fact none of them; 7 deals with non-collinear case which is not yours, and 5 inverts the spin, but this is again not what you want. You want one case with spin, and the other without. So it could be rather playing around with options 3 and 4. Option 4 smashes down the magnetism which is undesired, and option 3 allows to introduce spin moment on selected orbitals - one or more, in a row, and by repeating this you can generate any desired pattern of spins up or down. The use of DMtune is very straightforward - you start the script and answer its questions. However, you must be able to count your basis orbitals and know which one belongs to which atom and which (l,zeta-)states. Note that this will give you only the manipulated density matrix, which has no guarantee to survive in this spin configuration after some iterations. Good luck, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Chimie, Physique et Mat\'eriaux, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Question about spin configuration
On Sun, 30 Mar 2008, Cheng Pin Yen wrote: | Dear all : | I have one question about siesta. | I used siesta 2.0 version and system is a Ru | complex(Creutz-Taube ion) which has two redox site on | Ru atom. | One is 4d^6 : | t2g(dxy + dyz + dxz) was supposed to be 6 | electrons, | eg(dz2 + dx2-y2) was supposed to be zero . | And the other one is 4d^5 : | t2g was supposed to be 5 electrons , | eg was supposed to be zero . | How could I use command to do this ? Dear Cheng Pin Yen: as a rough tuning, you can define spin on each atom directly in the InitSpin block of fdf input. One of your configurations is magnetic and the other not, so with a bit of luck you will be able to stabilize one and another of them. As a more sophisticated tuning, you can hack the density matrix, addressing individual (on-site) terms for each orbital, and pray that these configurations survive in subsequent self-consistent calculation with a prudent mixing. I wrote a tool http://www.home.uni-osnabrueck.de/apostnik/Software/DMtune.tar.gz which, as it is, does not allow to modify occupations of individual orbitals, but it can be easily changed in this sense, I think. Independently on this technical problem, and speaking about Ru ion, I wonder how reasonable your results will be without a spin-orbit... Good luck, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Chimie, Physique et Mat\'eriaux, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] pseudo database;
On 13/03/2008, at 20:36, Andrei Postnikov wrote: | 1. Does somebody know whether there s a way to see from the .psf file if the | core correction is included, and with which radius? | On Fri, 14 Mar 2008, Eduardo Anglada wrote: | Just now there is no way, but I have all the info, of course you can ask for | it. Eduardo: would it be too much to ask you to (gradually) add this information to the comment line in each pseudopot ? :-) Best regards Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Chimie, Physique et Mat\'eriaux, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] DZP basis for elements with semicore states
On Wed, 23 Jan 2008, X. Feng wrote: | Hi Andrei, | | Thanks a lot for the explainations and the example for Cr. | | SIESTA manuel mentioned that to generate a polarization | orbital one needs to use P in PAO.Basis. Right | However, I saw | from Siesta website that in some of the optimized basis sets, | for example Ti and Mn, only a normal single zeta (4p) is used | as a polarization orbital. Then it is probably not a polarization orbital in a strict sence, derived from an s-orbital distorted by electric field. Rather, this is an orbital of p-angular symmetry, with radial dependence obtained or tuned I do not know how. | It is mentioned in the notes for these | bases that this normal single zeta is intended for polarization. What fo you mean by intended for polarization? Who intends to polarize whom? Can you give a reference? | Of course, this is not the definition of | a polarization orbital in the SIESTA manuel. No, it isn't. But from the point of view of Siesta, there are no fixed rules as for how to construct basis. You can draw you basis function by hand (, difitalize it) and go ahead with calculations. The question is, how good this basis would be. The standard scheme (DZP, Energy shift) is merely a reasonably good and fool-proof algorithm. | I don't know the reason | why people use a normal single zeta rather than following the manuel. Because they either believe that their tuned basis is more efficient, or they specifically want to use minimal basis, and they need to tune its functions as good as they can. | Is it better in some cases? Or, it is not a justified method? Depends on your definition of better and on cases in question. What do you want to achieve and at which price? The optimization of basis sets in any method using atom-type orbitals is a bit problematic. | By the way, do you have a tested basis sets for Cr, V? I did not tune them myself, if that's what you mean. However, I test basis sets (either standard, or borrowed from people) with respect to the tasks I have in mind for these elements. | I have pseudopotential for Cr and V without semicore states. | The magnetic moments are too high, so semicore states are necessary. This is a wrong argumentation. 3p semicore in Cr, V is not likely to affect your magnetic moment in a noticeable way, because 3p remains fully occupied. You'll see this if you compare total magnetic moments. So if your magnetic moments are too high it must be due to a different rason. Of course you may see that the figures for local magnetic moments, according to Mulliken populations, vary a bit as you vary basis functions, but this is due to ambigous definition of local magnetic moment, and has nothing to do with semicore as such. On the contrary, the presence of semicore may be quite important for comparing energies of different phases (e.g., magnetic phases), for getting right relaxation lattice constant (say, you may be 2-5% off the exp. volume with semicore and 15% off without it). | I'm currently testing PP's for Cr and V from a publication with | DZP, which will be generated by the method I learned today. | I'll be very grateful if you can send me PP's and basis, if you | have some. There is a PP + basis repository at the Siesta web site, and a lot of previous enquiries in the mailing list. Moreover, usually in their publications people describe what basis and PP they used. So if you want to profit from somebody's particular tuning for a particular system, you can write to people directly. You see, even for your Cr and V, the basis (and particularily whether you need semicore or not) depends so much on what you need (band sructure? phase diagram? phonons? transport?) and in which systems (bulk metals? clusters? organometallic molecules?) Good luck Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Problem FixSpin
Dear Marcos, the systems Eliseo is working on are molecular magnets, they usually have well localized spins on 3d ions and can be brought to self-consistency in a number of different magnetic configurations. Even if (usually) there is a certain (single-determinant) ground-state magnetic configuration, it is Eliseo's good right to expect a number of other metastable states to converge as well. The exact numbers (from integrated charge density or Mulliken analysis) are usually unimportant, because they reveal in one or other way a certain arrangement of integer spins of 3d ions. And of course the fixed spin moment, if imposed, remains fixed at each iteration (at the expense of introducing two Fermi levels, also in each iteration). Best regards Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+ On Wed, 28 Nov 2007, Marcos Verissimo Alves wrote: | Eliseo, | | The numbers you are showing for the spins are extracted from some | integration over the densities, or those are the mulliken charges on the | system? I will suppose that those are the numbers that come out after one | or more complete self-consistency cycles; if it is so, it means that your | system doesn't want to keep the initial values you have supplied; it | will be more advantageous, energetically, to have the final configuration | of spins. The initial configuration is just a starting guess, but | self-consistency takes over after this. Nothing guarantees that you will | have the final global net spin, or the individual net spins of the atoms, | equal to the initial values after a whole SCF cycle. Actually, you can | constrain the total magnetization of a system in siesta, but I am not sure | if the total magnetization would be constrained over all the scf steps. | | Cheers, | | Marcos | | | | Vous avez écrit / You have written / Lei ha scritto / Você escreveu... | Eliseo Ruiz | Dear Colleagues, | | I am calculating a periodic system with 40 Cu(II) cations in the | unit cell. I am trying to | calculate the energy for the different spin distributions. However, | only in a couple of | cases I obtain the signs of the spins that I specified in the input | file. For instance, for one | of the spin distributions: | | | # Spin defintion | SpinPolarized T | FixSpin T | TotalSpin 8.0 | %block DM.InitSpin | 1 +1.0 | 2 -1.0 | 3 -1.0 | 4 +1.0 | 5 +1.0 | 6 +1.0 | 7 -1.0 | 8 -1.0 | 9 +1.0 | 10+1.0 | 11+1.0 | | 40+1.0 | %endblock DM.InitSpin | | the results for the spin density: | atom up down total | 1 4.968 -5.401 -.433 | 2 4.978 -5.411 -.433 | 3 5.397 -4.948 .449 | 4 5.406 -4.942 .464 | 5 5.431 -4.963 .468 | 6 5.407 -4.962 .445 | 7 5.415 -4.977 .438 | 8 4.957 -5.391 -.434 | 9 4.943 -5.404 -.461 | 10 5.439 -4.964 .475 | 12 5.420 -4.987 .433 | | | The total spin is satisfied but not the local spins. I tried several | options of DM.MixingWeight and | DM.NumberPulay as well as for ElectronicTemperature but I cannot | reach the selected spin | distribution. | | best wishes, | | Eliseo | | | | -- | Dr. Marcos Verissimo Alves | Post-Doctoral Fellow | Unité de Physico-Chimie et de Physique des Matériaux (PCPM) | Université Catholique de Louvain | 1 Place Croix du Sud, B-1348 | Louvain-la-Neuve | Belgique | | -- | | Gort, Klaatu barada nikto. Klaatu barada nikto. Klaatu barada nikto. | |
Re: [SIESTA-L] Problem FixSpin
Dear Eliseo, with respect to your problem of initializing different spin configurations, I can offer two comments: 1. When you start from scratch, I am not sure how exactly the initial spin splitting is constructed; I also had sometimes an impression that the resulting splitting (in the 1st iteration) is too small. My rule of thumb in similar cases is to use maximal splitting, i.e., %block DM.InitSpin 1 + 2 - etc. in all cases, even for Cu. 2. As you have found a reasonable solution for one spin configuration you can invert certain spin(s) in its density matrice using my tool http://www.home.uni-osnabrueck.de/apostnik/download.html - DMtune that will normally give you a better starting point for further convergence than starting another spin configuration from scratch. Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+ On Wed, 28 Nov 2007, Eliseo Ruiz wrote: | Dear Colleagues, | | I am calculating a periodic system with 40 Cu(II) cations in the unit | cell. I am trying to calculate the energy for the different spin | distributions. However, only in a couple of cases I obtain the signs | of the spins that I specified in the input file. For instance, for one | of the spin distributions: | | | # Spin defintion | SpinPolarized T | FixSpin T | TotalSpin 8.0 | %block DM.InitSpin | 1 +1.0 | 2 -1.0 | 3 -1.0 | 4 +1.0 | 5 +1.0 | 6 +1.0 | 7 -1.0 | 8 -1.0 | 9 +1.0 | 10+1.0 | 11+1.0 | | 40+1.0 | %endblock DM.InitSpin | | the results for the spin density: | atom up down total | 1 4.968 -5.401 -.433 | 2 4.978 -5.411 -.433 | 3 5.397 -4.948 .449 | 4 5.406 -4.942 .464 | 5 5.431 -4.963 .468 | 6 5.407 -4.962 .445 | 7 5.415 -4.977 .438 | 8 4.957 -5.391 -.434 | 9 4.943 -5.404 -.461 | 10 5.439 -4.964 .475 | 12 5.420 -4.987 .433 | | | The total spin is satisfied but not the local spins. I tried several options | of DM.MixingWeight and | DM.NumberPulay as well as for ElectronicTemperature but I cannot reach the | selected spin | distribution. | | best wishes, | | Eliseo | |
Re: [SIESTA-L] question about energy of a single atoms
On Tue, 4 Sep 2007, bipul rakshit wrote: | hello siesta user, | i am doing calculations on Sn (tin)... | i do the convergence test for mesh cutoff and monhorsk pack. | the cubic Sn converges at 8x8x8 and mesh cutoff of 268 Ry. Since i use 8 atoms per unit cell so i get 512 atoms (auxillary supercell 4x4x4) | | Now my question is | 1. when i run siesta with the above parameters and at the end of the file the energy it give is the energy of 512 atoms or what Dear Bipul, the total energy is given per the number of atoms you really have in your cell, i.e. 8 in your case. The auxillary supercell is an artificial construction to facilitate keeping trace on overlaps across the cell boundary, in systems with periodicity. | 2. If its the energy of 512 atoms then how to find the energy of single atom, in the same volume occupied by a single atom?? To get total energy per atom, you divide your result by 8. I don't quite understand in the same volume occupied by a single atom. All all 8 atoms structurally equivalent? Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Fwd: [SIESTA-L] net charge calculation
On Thu, 19 Apr 2007, Yurko Natanzon wrote: | Dear Adam Gali, Andrei Postnikov, | thank you for explanation. Actually, I'm not interested in particular | numbers, but want to observe changes in ionic charges of atoms, | neighboring to dopant and the dependence of such changes on a dopant | concentration (comparabale to undoped system). Dear Yurko: yes, this is (the only) reasonable approach - to look at trends throughout a row of otherwise similar systems. | Is it possible to do it with Mulliken analysis? Yes (as with any other method to evaluate these charges) | And how does the charge depend on the basis? If you ask about systematic observations - I don't know. My general impression is - the more complex the basis, the more weird the numbers. I can only add that, should you be interested in getting charges integrated over atomic spheres, I've written a simple tool for this, which you can download from http://www.home.uni-osnabrueck.de/apostnik/Software/grdint.f You run it as a standalone code, using the RHO calculated by Siesta. The basis shape then plays no role, and you can compare the numbers with results from other codes, which use integration over spheres. (Of course the numbers do depend on the sphere radii). Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] net charge calculation
On Thu, 19 Apr 2007, Yurko Natanzon wrote: | Well, this works for TiO2, but when I tried to do it with SiO2, I've got: | | then oxygen has positive ionic charge (6-5.987), and silicon has | negative (4-4.026). Does it make any sense?? Dear Yurko: they have no more sense than ionic charge in general; such numbers are totally dependent on their definition (in the case given - on your choice of basis functions), even more so as you deal with covalent bond. Similarly funny numbers may occur e.g. in LMTO-ASA (charges come from integration over overlapping atomic spheres). It doesn't mean that your calculation is wrong, though. Just, better not to show these numbers to chemists... Best regards Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] force and energy optimization
Dear Saswata: your calculation is not converged. Reduce DM.MixingWeight; try DM.NumberPulay 2 Moreover check if you'd need more than one k-point, in order to get convergence of forces. Good luck, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] -- http://www.home.uni-osnabrueck.de/apostnik/ --+ On Wed, 28 Mar 2007, Saswata Bhattacharya wrote: | dear frnds, | i have got a problem where with the siesta run i have optimized force on the atoms very easily,but the total energy is not optimized.this should not be the normal case as force is the derivative of energy..i am attaching both the input .fdf file and output .out file so that with these u can make comment and with this i will be greatly benefitted.please tell me what do i need to do to get the energy optimized. | regards, | Saswata | | | - | Heres a new way to find what you're looking for - Yahoo! Answers
Re: [SIESTA-L] DOS
On Tue, 6 Feb 2007, Neil Dixon wrote: | dear all, | i am a new user of siesta.can u please tell me which command give me systemlabel.DOS file in my output.I have set the range of PDOS and LDOS but yet i have not got any systemlabel.PDOS ya systemlabel.LDOS files.is there any command which gives the respective files in the outpit? | please hel | with regards, | Neil Dear Neil: PDOS and LDOS are different properties which are conrolled by two different blocks. In order to produce .LDOS as r-resolved function on the grid, use e.g. %block LocalDensityOfStates -2.30 -1.30 eV %endblock LocalDensityOfStates in order to get .DOS and .PDOS (functions resolved in energy), use e.g. %block ProjectedDensityOfStates -27.0 13.0 0.2 800 eV %endblock ProjectedDensityOfStates I don't think you need any other switches to activate either one, or another. What may be a reason for missing files in you case: the calculation of LDOS / PDOS is only done after the electronic iterations loop is normally terminated, therefore you have to achieve the convergence, or limit the number of iterations. Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Question about LDA and GGA pseudopotential
On Fri, 2 Feb 2007, Dr. Y. Lin wrote: | BTW: What is the MeshCutoff? I read the manual, and it says: | | MeshCutoff | (real energy): Defines the equivalent plane wave cutoff for the grid. | I thought siesta is not using plane wave basis. Why does the code needs this? | Could somebody point out where I shall get more info? Dear Jin Zhang: chack the basic Siesta paper JPCM 14, 2745, Sec.6 (The paper is algo good for other questions) Concerning the previous discussion, LDA vs. GGA bond lenth, the effect of grid etc: the fine real-space grid is needed in order to stabilize, i.e., to converge (in terms of grid cutoff) your results like total energy, forces and everything that is related to them, including bond lengths. However, a good grid offers no fix for everything else and would not repair what was due to bad pseudopotential, or basis, or both. For system like CH4 one can easily make a very accurate test for LDA truth or GGA truth with a good all-electron code (or, such benchmarks are certainly available), and then you'll see whether your deviation from experiment is a known LDA fault, or a fault in your SIESTA setup which fails to reproduce the LDA truth. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ | | Thanks. | | | Jin Zhang | | On 2/2/07, Dr. Y. Lin [EMAIL PROTECTED] wrote: | | On Thu, 1 Feb 2007, Oleksandr Voznyy wrote: | |1) Why is GGA pseudo-potential gives worse result than LDA? |What do you mean worse? |A well converged calculation (i.e. converged in E cutoff, k-grid, and |converged vs pseudo cutoff radii as well) |LDA SHOULD give bonds (or lattice constants) 1% smaller than expt (while |binding energies are strongly overestimated) |while GGA usually overestimates bonds by about 2% (while binding | energies |(for molecules) are much closer to expt). | | First, thanks for the info. Let me explain why I think GGA gives worse | result. | | Let's take CH4 as an example. Here's the fdf input file I have | - Start of the file --- | SystemName Methane Broyden Optimization | SystemLabel general | NumberOfAtoms 5 | NumberOfSpecies 2 | | MeshCutoff 50 Ry | | %block ChemicalSpeciesLabel | 1 6 C # Species index, atomic number, species label | 2 1 H | %endblock ChemicalSpeciesLabel | | #PAO.BasisSize SZ | %blockPAO.BasisSizes |C DZP |H DZP | %endblock PAO.BasisSizes | | AtomicCoordinatesFormat Ang | %block AtomicCoordinatesAndAtomicSpecies |1.0582 0.9353 0.8103 1 |1.4145 1.5662 0. 2 |1.2065 1.4452 1.7588 2 |0. 0.7294 0.6710 2 |1.6121 0. 0.8114 2 | %endblock AtomicCoordinatesAndAtomicSpecies | | XC.functional GGA | XC.authors PBE | PAO.EnergyShift 0.02 Ry | | DM.Number.Pulay 3 | | WriteForces T | | MD.TypeOfRun Broyden | MD.NumCGSteps 400 | MD.Broyden.History.Steps 6 | MD.Broyden.Initial.Inverse.Jacobian 1.0 | --- End of the file --- | | This input file is used along with GGA pseudo potential of Carbon. | For LDA pseudo potential case, I just change the line with Carbon | at block PAO.BasisSizes to | | CDZ | | Basically, this setting reduces the total basis of Carbon atom from 13 to | 8. | | Then I ran siesta to relax the CH4 and measure C-H bond length after | relaxation: | | Exp. LDACarbon GGACarbonCCCBDB calculated |PBEPEB/6-331G** | 1.094 1.097 1.10 1.0977 | | It seems LDA pseudo potential did give a shorter bond-length while still | longer than experimental value. GGA pseudo potential gives an even longer | bond-length, which is obviously, a worse result. | | |2) Is it safe to use LDA pseudo-potential while setting GGA-PBE as |x-correlation in input.fdf? |No, it is not recommended. |However, you can use the cutoff radii from LDA pseudo to generate the | new GGA |pseudo. Have you tried to compare 2 input files that you have? | | | I guess it is not very good to use LDA pseudo-potential while setting | GGA-PBE as x-correlation. Thanks for confirming this. However, I'm not | quite sure I understand your suggestion. | | Let me clarify this: I'm not trying to generate pseudo-potential myself. | So, all the input I said are for Molecular Dynamics or Relaxation. Under | this situation, I'd say, basically, my input files are the same except for | the basis sets. For GGApseudo, I set
Re: [SIESTA-L] HOMO
On Mon, 15 Jan 2007, navaratnarajah kuganathan wrote: | Dear Users, | | I would be thankful to you, how can i generate HOMO and LUMO orbitals for the C60-porphyrin complexes(nearly 600 atoms).Thanks in advance Dear Kugnathan: Calculate Wavefunctions using e.g. %block WaveFuncKPoints 0.00 0.00 0.00 from your HOMO to your LUMO %endblock WaveFuncKPoints ( make sure that you have SolutionMethod T ) and then vizualize them with Denchar. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Fermi level and k-grid
On Wed, 10 Jan 2007, Oleksandr Voznyy wrote: | The Fermi level is normally calculated by setting the cumulative occupation | number of all bands to the number of valence electrons. | | As I understand this means that Ef in semiconductor would always be | at the VBM and not in the middle of the gap? Alexander - sorry, it seems that I was wrong. I just checked my old calculation for a wide-gap dielectric and see that Fermi energy = -6.866874 eV while the energies bordeing the gap -8.08 and -4.66. One should look into details of implementation... Yes in other band structure codes I know the Fermi energy is fixed by the last occupied band, i.e. it is set at the valence band top. Then if it technically lies higher, this must be due to the energy broadening introduced. But apparently in Siesta it is done differently. | By the way, how the bandstructure is calculated using only several k-points? The band structure as such (continuous bands) is not calculated, each k-point enters independently of others and contributes a (smeared) peak in the density of states. This summary DOS is (roughly speaking) integrated, and as the number of electrons is achieved, the Fermi level is set. (The details of implementation might be different). Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Fermi level and k-grid
On Wed, 10 Jan 2007, Oleksandr Voznyy wrote: | Hi, | till recently I though that checking the convergence of total energy vs k-grid | cutoff is enough. | However, now I've found that while total energy can be very well converged, | Fermi level position is not, and requires at least twice denser k-grid (and ~4 | times more time). | | Here is my example for bulk GaAs: Dear Alexander, the problem you describe stems from the fact that k-space summation, done by sampling, is not accurate enough. However, this is technically no problem for systems with large enough band gap and/or molecules. On the contrary, for metallic systems, or those where the Fermi level crosses states in the gap, the convergency of results with k-points in indeed disappointingly slow (as compared with codes using tetrahedron integration). The Fermi level is normally calculated by setting the cumulative occupation number of all bands to the number of valence electrons. As I understand the situation in Siesta, this cumulative occupation is obtained by integrating the Fermi function smeared with ElectronicTemperature, and summed up over k-points with their respective weights. That's why increasing ElectronicTemperature usually suppresses the fluctuations and helps the convergency, but somehow deteriorates the resulting energy values. Considering your example, for pure GaAs (or any semiconductor) you probably won't care much about the exact value of E-Fermi because you get the valence and conduction bands right, and total energy is stable. Etot is more stable because it is integral property while the calculated E-Fermi is differential one which shifts back and forth with every single k-point added. Correspondingly, the DOS is a sum of smeared delta-peaks, it also wildly changes with adding k-points, and converges extremely slowly to the DOS found from other band structure code with tetrahedron integration. In a metal system, metal Etot won't be so stable against k-mesh as in semiconductor, and would normally require much much more dense k-mesh. It would be such a good investment to add a tetrahedron integration in Siesta... And a not very difficult one... Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ On Wed, 10 Jan 2007, Oleksandr Voznyy wrote: | Hi, | till recently I though that checking the convergence of total energy vs k-grid | cutoff is enough. | However, now I've found that while total energy can be very well converged, | Fermi level position is not, and requires at least twice denser k-grid (and ~4 | times more time). | | Here is my example for bulk GaAs: | kgridEf, eV k-pnts SCFtime forces Etot | cutoff | 8 -5,3105 32 1 0,00509 -789,50149 | 10-4,9698 -- -- 2,19E-4 -789,44567 | 12,2 -4,1639 108 -- 0,0052 -789,50521 | 16,287-4,1839 256 1,370,0028 -789,50544 | 20,359-4,7579 500 -- 1E-3-- | 26,467-5,1201 11833,139 2E-4-789,5052 | 30,5 -4,9783 1800-- 4E-5-789,5054 | 32,575-4,7802 20484,736E-6-789,50545 | 40,7 -4,7676 4000-- 1,28E-4 50 -5,1057 16 | 2,07E-4 | | As you can see, Fermi level varies in the range of 0.6 eV!!! (all the bands | don't shift). The middle of the gap is at -4,75 | | My questions are: | | 1. How Fermi level is calculated? | Is it just filling the available bands with a given amount of electrons (based | on calculated DOS on a given k-grid) after all calculations are done? What | smoothing of DOS is used then??? | | Ef doesn't converge actually. | If 1. is true then I can understand it - DOS shape changes quite significantly | and real convergence would be only when one gets all possible k-points. | | 2. Since the total energy is calculated on the same k-grid, why it doesn't | show the same behavior? i.e why it is less sensitive than Ef? | Should one bother at all about Ef during geometry relaxation? | | 3. Does the filling of the bands affect the forces on atoms, and thus | explicitly affects total energy? | Imagine such a situation: | a molecule adsorbs to a dangling bond on a surface only if it is empty or only | partially filled, | if I set the Ef 0.5eV higher, I make the dangling bond completely filled and | moleculed would not adsorb at all, | i.e. we end up in a completely different geometry and total energy. | | I will appreciate very much any comments or suggestions. | | Sincerely, | Alexander. | |
[SIESTA-L] DM manipulation tool
Dear Siesta community, I offer to you attention an updated tool (attached tar) which allows to manipulate (on a very elementary level) the density matrix, in order to perform several tasks which arise from time to time, e.g. C Possible use: C 1) read DM unformatted, write formatted, or C 2) read DM formatted, write unformatted C-e.g., for porting to a different system; C 3) expand non-magnetic (nspin=1) to magnetic (nspin=2) case, Cintroducing magneting splitting on selected atoms; C 4) contract a magnetic case (nspin=2) into a non-magnetic one (nspin=1) C- e.g., for more economic continuation of a case Cwhich has lost magnetism in the course of calculation; C 5) flip spin on a selected atom, to make a (hopefully) Cgood start for a different magnetic configuration. C 6) expand collinear DM (nspin=2) to non-collinear format (nspin=4), Cand/or to include spin-orbit (nspin=8), C 7) rotate non-collinear DM on selected atom. Many options have been recently added, and may contain bugs. Please feel free to use and report to me any problems. No documentation is provided, because one needs simply make and start dmtune - the rest is (hopefully) self-explaining. Please be advised that the version with spin-orbit (nspin=8) is not supported by siesta 2.0, and for the details on this version you should address J. Ferrer and L. Fernández-Seivane, as has been mentioned several times in the list. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ DMtune.tar Description: Unix tar archive
Re: [SIESTA-L] Why it is so hard to converge with 32 atom unit cell.
On Tue, 5 Dec 2006, Vasilii Artyukhov wrote: | Finally, a word regarding the k-point grids: although a dense k-point grid | is essential for an adequate description of properties of metals, I don't | think that denser grids could help you improve your SCF convergence, since | the SCF equations are solved at each k-point independently (which is why the | k-point parallelization is so efficient). A comment on this: in my understanding, not the SCF equations are solved at each k-point independently, but the Kohn-Sham equations which are only part of the SCF loop (if we discuss the conventional diagonalization scheme). Then, in order to find the Fermi energy and re-define the density matrix, the information is needed over all k-points. If the k-mesh is sparse - and the system is metal - then a too high weight attributed to a given k-point, in which a given band might happen to be either occupied, or empty, would lead to additional fluctuations/instability. More dense k-mesh reduces the importance of every single k-point and thus stabilizes the convergence. Of course by intruducung a large enough artificial smearing (Electronic Temperature) one can suppress this instability, at the expence of overall precision. Probably that's OK in order to enforce some initial convergence. The exact behavior and the balance of these factors is of course up to a system in question. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Why it is so hard to converge with 32 atom unit cell.
+-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ On Tue, 5 Dec 2006, siesta88 wrote: | Dear SIESTA users: | I used dzp atomic basis,160 ry meshcutoff ,30 Bohr kgrid cutoff to calculate FePt system(see the pictures of my unit cell in attachment), however, it seems to take forever to converge with the SCF cycle. dDMax(around 1.0) stays far beyond tolerence (0.1), and if I increase the cell to 108 atom unit cell(3x3x3 supercell), dDmax becomes even larger. Yet the simplest unit cell with 4 atoms converges fine. | | Here is a fragment of my SCF cycle: | siesta: 42 -25183.5213 -24873.4575 -24873.4883 1.2337 -2.1597 | siesta: 43 -25153.3555 -24872.5874 -24872.6213 1.0953 -1.7754 | siesta: 44 -25204.6434 -24870.7750 -24870.8201 1.0794 -2.1741 | siesta: 45 -25174.8238 -24868.9171 -24868.9419 1.3392 -2.3371 | siesta: 46 -25138.7377 -24868.4595 -24868.4879 1.0200 -2.1080 | siesta: 47 -25177.5378 -24866.3887 -24866.4199 1.1245 -1.9894 | siesta: 48 -25159.6251 -24865.8291 -24865.8733 1.0482 -2.0247 | siesta: 49 -25206.9572 -24864.1873 -24864.2264 1.0436 -2.1378 Dear SuiYang: Your mixing parameter is too large redata: New DM Mixing Weight = 0.2500 that's why (at least) you get no convergence. Try to set DM.MixingWeight.05 (or less), to begin with. There are two important issues to consider: 1) in a large supercell, there is genererally more difficult to achieve convergence because you have many bands with contributions from similar atoms around the Fermi level. So the system has a tendency to swap charge from one atom to another, which effect you might need to damp more eficiently than was the case in a smaller cell. In the single cell with only few bands and few atoms, you won't normally have this degree of freedom in your system, so the convergency goes sraightforwardly. 2) In a metal system, depending on the actual band sructure, the accuracy of the k-summation/interation becomes an important issue influencing the accuracy of your result, and the convergence to it. In SIESTA where no tetrahedra integration is available, you must carefully test the effect of increasing number of k-points, and the energy broadening used in the search for the Fermi level. Your present choice 5*5*5 k-points with ElectronicTemperature 300 K is not a priori unreasonable for getting a preliminary convergence, but be careful to make additional tests if you want to extract some energy differences afterwards. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] WrieWaveFunctions
On Wed, 29 Nov 2006, Cherry Y. Yates wrote: | I was using siesta to get wavefunctions. From the | manu, if WriteWaveFunctions is .false., it will | generate a .WFS file, which can be read by the utility | READWFS, but I can not find this file. If | WriteWaveFunctions is .true., it will write | eigenvectors in the main output file, but I cannot | find this part in the output file either. Please let | me know how to get this infomation. Dear Cherry, you also need to define WaveFuncKPoints - by default this list is empty, that's why you get no output. | I also wonder if | siesta provides utility to draw molecular orbitals, | HOMO LUMO, etc. Yes, with Denchar. See the documentation in Util/Denchar and make executable in the main Siesta source directory with make denchar. Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Phonons in CdSe
On Tue, 21 Nov 2006, marcel wrote: | Dear All, | | I'm doing some test calculations on bulk hexagonal CdSe. I want to calculate | the phonon dispersion. | (Im using the version of T. Archer, for LO-TO splitting and include BC in my | Siesta calculation. | The dielectric constants I use are from literature on measurements. ) | | I get perfect acoustic frequencies at 0.0 cm^(-1). | Lowest optical and highest optical frequency are underestimated by max. 2 % | from experiment. OK. | However, the mid-range frequencies are overestimated, up to 10%. | | So would you think this is ok, or should I further investigate the source of | this discrepancy? Marcel: it can be deduced from your message that you have this problem already in Gamma, and for the minimal cell. Did you check the effect of k-mesh? The forces and force constants may vary with it. Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] atom label ordering problem
On Tue, 21 Nov 2006, Chu Chun FU wrote: | Dear Postnikov, | | I could but, I didn't since in principle the result should not depend on the | atom label ordering... | | Chu Chun Dear Chu Chun, What I mean is, if you for instance always set initial spin for type 2 only: if your type 2 is Fe it converges to magnetic solution as it should be, and if your type 2 is Cu it doesn't manage to induce enough spin on Fe and converges to a non-magnetic solution. Andrei
Re: [SIESTA-L] atom label ordering problem
On Tue, 21 Nov 2006, Chu Chun FU wrote: | Dear All, | | We are studying the FeCu system with Siesta and we have found some problems | related to the | labels (1 or 2) asigned to each chemical species. | | For example when considering the system containing 1 Fe atom and N-1 atoms of | Cu with BCC | structure (e.g. N = 54), the results (total energy and magnetic moment of Fe) | differ depending | on which atom is labeled as atom '1': | | If we put | | %block ChemicalSpeciesLabel | 1 29 Cu | 2 26 Fe | %endblock ChemicalSpeciesLabel | | we obtain reasonable values for total energy and magnetic moment for the Fe | atom (about 3 | mu_B). | | However when we change the order and put | | %block ChemicalSpeciesLabel | 1 26 Fe | 2 29 Cu | %endblock ChemicalSpeciesLabel | | we obtain a different value of total energy (0.6 eV higher) and zero magnetic | moment for Fe ! | | We wonder if anybody had a similar problem before ? or any idea on the origing | of this problem | ? Dear Chu Chun: and do you set %block DM.InitSpin differently, depending on which type is Fe ? Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] About Atomic Forces.
On Tue, 21 Nov 2006, SuiYang wrote: | Dear SIESTA users: | I have used SIESTA to study the property of FePt systems.Here I have a question about the atomic forces SIESTA calculated: | When I am calculating the FePt unit cell(including 4 atoms) with its lattice constant around 4 ang.The corresponding atomic forces are around 0.02 eV/ang.When I am using the 2x2x2 supercell(including 32 atoms) as unit cell, with relatively the same atomic positions,and the lattice constant is changed to 8 ang to maintain it is the same bulk as the previous one,however ,the atomic forces are increased to around one hundred eV/ang. | | I also calculated another system using 3x3x3 supercell as unit cell, the forces are much even higher ,it seems the atomic force will increase with the number of atoms and the lattice constant. | | I am rather confused with its behavior,why there is such a drastic change in atomic forces and stress tensor while the bulk is the same, only the expression of it is different. | | I've attached both 1x1x1 and 2x2x2's input and output files,if needed. | I really need your help. | Dear SuiYang, As Riccardo pointed out, you are doing calculation with the default setting for k-points, that is just one k-point. This is no good when calculating bulk - in any case, metals. You can see that Siesta fails to converge your system - which is in fact metal with partially filled bands, - when using just one k-point and hence no dispersion; it just desperately swaps the Fermi energy back and forth: siesta: iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV) ... siesta: 23-3027.3747-2953.1394-2953.2110 0.8298 2.5907 ... siesta: 98-3034.0919-2952.4118-2952.4146 0.8258 -2.0656 siesta: 99-3027.3736-2953.3750-2953.4466 0.8258 2.5911 Normally you should be wondered by this behavior before going any further... After you manage to get reasonable convergence of electronic properties, take into account that the forces are especially sensitive to the k-mesh, and (if you really need good forces, i.e. for phonons) you should check how they converge as you increase the number of k-points. Your forces are zero in your 4-atom cell, because the positions of atoms are symmetric with respect to the real-space mesh (defined by the Mesh cutoff, 200 Ry in your case), and the errors cancel. In order to see how good/bad the forces REALLY are, you can displace the atoms a bit from the symmetric positions, imposing the uniform shift (AtomicCoordinatesOrigin). In your 2x2x2-cell you have apparently the atoms situated at random with respect to the mesh along the Y-direction, and you see huge fluctuations of forces, which - in principle, i.e. apart from numerical errors - should be zero by symmetry: siesta: Atomic forces (eV/Ang): 10.01 94.7209640.06 2 -0.04 -94.7209750.03 30.04 -117.510866 -0.02 4 -0.05 117.5108590.00 I must confess I don't understand why X and Y behave differently, avan as they are completely equivalent with respect of atomic positions and mesh construction. I'd be glad to see somebody's comment on that. You did't need to provide .fdf separately because they are dumped at the beginning of each Siesta output anyway. I don't know what you are after in your study, but if I remember correctly, in real life FePt is magnetic, so that when you switch on the spin, this will affect your forces as well. Hope this helps, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] HCP lattice vectors
Emma, the suggestions done by Alexander on your problems seem very reasonable. I'd address some more points. 1. In the pseudopot generation, you use Co 3s2 3p6 3d7 i.e. for the Co2+ case, and the basis functions 3s(SZ) 3p(SZ) 3d(DZ) 4s(DZ) 4p(SZ) I wonder if you really need 3s, 3p as valence state. You may wish to include them in a FLAPW code (along with I4d) when this would cost you nothing, but in a pseudopotential calculation you have to decide which are your valence states needed for the chemical bonding, and I think 4s is much better bet. This will give you a much less problematic pseudopotential. In fact Co3s are too deep to have an effect on anything. Co 3p are not so critical, you can include them in the valence band if you wish, but this won't be really needed, because they are much lower than I5s. I think you could make better pseudopot with, say, 4s1 3d8 configuration; with core correction. For iodine, the pseudopot configuration seems reasonable, but I don't know about the details, radii etc. In the basis, I don't think you need 6s and 5d, better to add polarization orbitals on 5s and/or 5p. 2. Apart from this, there must be a minimum of total energy vs. volume, if you have good overlap of basis functions. But this minimum might be quite apart from the experiment, due to problems with pseudopots/bases. In your system, I'd expect the compressibility to be very different in the XY-plane and with respect to Z-compression. Try to explore the energy minimum within a large enough interval of varying a, and then c, or just start an automatic relaxation as Alexander suggests. 3. For the true relaxation, you'd need low force sum, that means higher EnergyCutoff. This will however depend on the choice of pseudopotentials. Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ On Thu, 9 Nov 2006, Emma Sceats wrote: | Andrei | | Thanks for taking the time to respond to my email. | | I must apologise for a couple of typos in my email to you. First off, I meant | to say that CoI2 crystallizes in P-3m1 (space group 164) rather than P3-m1. My | coordinates should have read: | | Co0 0 0 | I 1/3 2/3 0.25 | I 2/3 1/3 0.75 | | (oopps - big typo there!) | | I specify the lattice constant (3.96 Ang) and the Lattice parameters as: | | 1.0 1.0 1.679 90 90 120 | | I have attached the pseudopotentials and output file of a calculation. I | should say that I have successfully used the Co pseudo and basis set in | calculations on CoO (I chose this as a test system since there is lots of data | (both experimental and calculated) on the structure, bulk modulus etc and it | is quite similar (in terms of level of covalency) to the CoI2 system I am | interested in). Also the I pseudo and basis were used in some calculations I | ran on KI (and gave a good lattice const and bulk modulus). | | A couple of things to point out: | (1) the forces on the atoms in the unit cell are HIGH (0.9 eV/Ang) but only in | Z. | | (2) when I vary the lattice parameter slightly (as if I was calculating the | cohesive curve to get the bulk modulus) I do NOT see the expected parabolic | dependence of the energy on the lattice constant (with a minimum at the | equilibrium lattice constant). Instead I see an increase in the energy with | increased lattice constant. (although I do see an approximately linear | decrease in the pressure with increasing lattice constant). | | I am completely confused as to why I see these effects, and emailed the | mailing list originally because I thought I had done something wrong with my | input coordinates. | | Can you check the coordinates and lattice parameters are consistent (I was a | little confused with what you thought I should change) and if they are OK have | you any other suggestions you can offer? (changing the | pseudos/basis/coordinates/lattice vectors or parameters etc etc) Anything you | think I could try would be helpful - I am runing out of ideas!!! | | | Emma |
Re: [SIESTA-L] Pseudopotentials for F
+-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ On Thu, 9 Nov 2006, Ö£ÈÊ»Û wrote: | Dear Prof. Andrei Postnikov: | nbsp;nbsp; when I pg.sh the following Fluorine Pseudopotential | with siesta-2.0, no vps or psf filenbsp;isnbsp;created and the errors | are:nbsp;cp: stat 'VPSOUT' failed; cp:nbsp;'VPSFMT' failed.nbsp; | Is there any methodnbsp;for solving it? -- Yes, check the format of pseudopot input file; it might be that something got displaced when I pasted the file into the message. Here the file again, in the attachment. Good luck, Andrei Postnikov f01.tm2.inp Description: chemical/gamess-input
Re: [SIESTA-L] MD-calculation
On Sat, 21 Oct 2006, Michael Shin wrote: | Hello Dear users | I am trying to find the Eqlubrium volume of my system. I found it from total enery calculations . now I want to find it again using MD-CG methohs. | During calculations the Pressure is fluctutating and is given below. | The volume seems to be constant. | Can any one tell me that is there is some problem in my calculations? Michael: most probably you forgot to set MD.VariableCell T Moreover (although not related to you problem) check that your MeshCutoff (150 Ry) is sufficient for your system (look that the sum of forces over all atoms be small enough), because it is essential for geting forces and relaxation right. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] probleme_with_phenol_input
| Error in Cholesky factorisation in rdiag | Stopping Program from Node:0 | | I don't know what is exactly the problem ? A good chance that it is due to an error in input, e.g. two atoms at the same place Good luck, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] fermi level
On Tue, 26 Sep 2006, Marcos Verissimo Alves wrote:
| Hi Andrei,
|
| Your answer to Swaroop got me relieved, because I wasn't getting the same
| for my spin-polarized calculations... However, regarding his question:
| then the Fermi level for a plot of the total DOS (DOS_{up}+DOS_{down})
| would be just the highest of the two Fermi levels, correct?
|
| Marcos
Hi Marcos:
of course the way of graphical expresion is very much
to one's tastes. However to my opinion the essence of playing around
with FixedSpin is exactly to maintain and underline the difference
betwen two spin channels, including the difference in their
respective Fermi energies. This latter difference can be associated
to the splitting in an effective magnetic field.
To this end, I don't immediately see why showing THE SUM of
two spin-split DOS, coresponding to the FSM case, would make
any sense at all... Could you please explain?
If you bear in mind a comparison to some specroscopy experiment,
then we should look into how the spectroscopic events in two
spin channels sum up in the total spectrum - and bring them (the events)
to the same energy scale. For X-ray emission this will be, say,
(3d-up states up to Fermi-up) - (2p-up) +
(3d-dn states up to Fermi-dn) - (2p-dn) ,
and you should take into account the magnetic splitting of 2p core states.
In the XPS, the binding energy in each spin channel will be measured
from its relative E_F, etc.
A special case is when you have a system with band gap
in each spin channel, and so that the gaps overlap. In this case you
can probably draw a common Fermi level through this common band gap
quite meaningfully, even if the HOMO energies are different...
Does this correspond to your situation?
Best regards,
Andrei
+-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+
| Paul Verlaine University - Institute de Physique Electronique et Chimie, |
| Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France |
+-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/
--+
Re: [SIESTA-L] fermi level
On Mon, 25 Sep 2006, Sairam Swaroop wrote: | Dear siesta users | | I am running a spin polarized calculation and get two fermi energies for | fermi up and fermi down. How do i calculate the fermi level for a DOS | plot. Is it the average of these two levels.. | | thanks | | swaroop Dear Sairam Swaroop: two Fermi energies appear in the output if you perform calculations with fixed spin moment FixSpin T (because each of spin-up and spin-down values remains constant in this case, not just their sum). Check if this is what you really want. If yes, then you'll probably want for each spin channel to show its corresponding Fermi level: the splitting of two Fermi levels kind of represents the artificial magnetic field which is imposed in order to enforce the fixed moment. The average of both E_F doesn't make ay sense. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Vibra problem
On Tue, 18 Jul 2006, bipul rakshit wrote: | hello siesta user, | | i am facing a problem in using Vibra for PtN calculations. | i got frequencies at gamma point 0.0032 0.0488 0.0540 526.0463 526.0609 526.0764 | | it is correct? Three acoustic frequencies are close enough to zero; that's OK. Three others are almost equal, that is also a good sign. As for the value you obtain (how good it is), check the literature... | but when i construct a supercell = 1 in all the three direction. and when i run the siesta command, then it says fail to converge, eigen value problem | | what can i do??? plz help me with this regard... Try so set Diag.DivideAndConquer to F as has been discussed many times in the list Hope this helps, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Confusion with output
On Thu, 13 Jul 2006, lan haiping wrote: | Hi, All! | | I have a confusion with siesta's output . | I have set the AtomicCoordinatesFormat to* Ang*, | But i found the atomic coordinates in output with the same values in *Bohr . | * | ** | Part of output is below, would you please help me ? Yes. You misprinted the label AtomicCoordinatesFormat for AtomiCoordinatesFormat Ang so your Ang entry was ignored, and Bohr taken by default. Best regards, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ | regards, | | Hai-Ping | | *Input :* | ** | *AtomiCoordinatesFormat Ang | %block AtomicCoordinatesAndAtomicSpecies | 0.000 0.000 0.000 1 | 1.888 1.888 4.743 1 | 0.000 1.888 2.372 1 | 1.888 0.000 7.115 1 | 0.000 0.000 1.973 2 | 1.888 1.888 6.716 2 | 0.000 1.888 4.345 2 | 1.888 0.000 9.088 2 | 1.888 0.000 5.141 2 | 0.000 1.888 0.398 2 | 1.888 1.888 2.770 2 | 0.000 0.000 7.513 2 | %endblock AtomicCoordinatesAndAtomicSpecies | * | *Output :* | ** | *siesta: Atomic coordinates (Bohr) and species | siesta: 0.0 0.0 0.0 11 | siesta: 1.88800 1.88800 4.74300 12 | siesta: 0.0 1.88800 2.37200 13 | siesta: 1.88800 0.0 7.11500 14 | siesta: 0.0 0.0 1.97300 25 | siesta: 1.88800 1.88800 6.71600 26 | siesta: 0.0 1.88800 4.34500 27 | siesta: 1.88800 0.0 9.08800 28 | siesta: 1.88800 0.0 5.14100 29 | siesta: 0.0 1.88800 0.39800 2 10 | siesta: 1.88800 1.88800 2.77000 2 11 | siesta: 0.0 0.0 7.51300 2 12 | | initatomlists: Number of atoms, orbitals, and projectors: 12 180 208 | | | *
Re: [SIESTA-L] Vibra in siesta 2.0
On Thu, 13 Jul 2006, bipul rakshit wrote: | hello siesta user, | i am doing the phonon calculation with the help of Vibra from the siesta2.0 package. | i am doing PtN with LDA approximation. But in the dispersion plot only two modes are seen. Why not all the modes are visible. | i am sending the input files used in the calculation | ptn.fdf for fcbuild and vibra | and ptn-siesta .fdf for siesta | can any body have a look and tell me where i commit mistake Bipul: as you give no output, I can only guess from your vibra input that you only get Gamma phonon because you don't define SuperCell and hence you won't get any phonon dispersion. You define many q-values for the plot, but they will probably give you always the same numbers for the frequency, because you don't have any information in your calculation about how the force constants change in the real space from neighbor to neighbor, and this is crucial to get the phonon dispersion, i.e. to Fourier-transform this real-space dependence of the force constants. Still, you can already check from your calculation whether your Gamma phonon is OK - i.e., you must get three nearly zero numbers and three nearly-equal numbers, which would the the triply degenerate Gamma-TO phonon in PtN; you can certainly find a reference to compare. If this is OK, and you need to calculate phonon dispersion, you'd need to construct a supercell large enough, i.e. SuperCell_1,2,3 0 to assure that the force constants fall down within this supercell. You ask about instructions to vibra - the documentation file is in Util/Vibra/Docs Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] DOs of molecules
On Tue, 4 Jul 2006, Michael Shin wrote: | Now I used the |%block DM.InitSpin |1 0.000 | %endblock DM.InitSpin. | I got 0.00 moment. | | It means Spinpolarization calculations is not a Ferromagnetic calculation? Michael - it's about diferent things. You'd switch spin polarization if you want to treat spin-up and spin-down blocks differently, i.e. you allow (and/or expect) different occupation of spin-up and spin-down in the density matrix and other properties. So this doubles your calculation time (always), irrespectively of whether your spin-up and -down properties are in fact diferent or not. Ferromagnetic etc. is about the orientation of local magnetic moments at different atoms. May be e.g. antiferromagnetic or other. | but for Test purpose I also used Fe where I didnt use | the %block DM.InitSpin and I just used Spinpolarization=T | and at the end I got the magnetic moment, | is it becasue that Fe is magnetic in nature. You got magnetic moment because without the DM.InitSpin block ALL atoms (i.e. one in this case) were initialized spin-up, that happens to be exactly what you want for ferromagnetic Fe. | If I have 100 C-atom in the unit cell so I will have to specify | 0.00 initial spin for every atom? No, read the Manual on the DM.InitSpin. Keep the block, but list there only the atoms which have to be initialized as spin-polarized. If you don't have such atoms probably you don't need spin polarization at all. For the future: please choose more specific Subject line when posting to the maillist (to facilitate subsequent browsing), because what we are discussing has nothing to do with DOS in molecules. Sorry. Good luck, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] DOs of molecules
Dear Michael, you got a high-spin configuration of the carbon atom: 2s electrons are spin-compensated but two 2p-electrons have parallel spins. It may be counter-intuitive but physically plausible, and fully in accordance with the Siesta input logic. You declared spin-polarized calculation but did not specify a block DM.InitSpin In this case (see UsersManual on DM.InitSpin) DM.InitSpinAF takes control by default and sets maximal-spin ferromagnetic configuration of ALL ATOMS for the start. The remedy depends on what you want to calculate. Probably you just want to have a trial spin somewhere in your system (e.g., 3d atom or other defect). Then use the block DM.InitSpin to initialize exactly whatever you want; the other atoms (not included in the block) will be non-magnetic on start. Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ On Tue, 4 Jul 2006, Michael Shin wrote: | Dear Siesta users | I am doing some calculations on Carbon atom using GGA as well as LDA. For Basis I used DZP and MeshCutoff 150. Ry. | I used spinpolarized calculations and at the end I got some strange result i.e its shows that Total Spin1= 3 and Total Spin2=1 | Is it means that C-atom has magnetic moment or there is some problem in my calculations. For LDA I used the same C.psf which was provided with Siesta but the magnetic moment is same in LDA and GGA. I also got DOS and see the attached DOS file. | Here is the part of my output. | I will be thnx to everyone if some one guide me because I have no experience with magnetic structure | | Thank you very much | Mic. | | SietaOutputGGA. | GHOST: No ghost state for L = 0 | GHOST: No ghost state for L = 1 | GHOST: No ghost state for L = 2 | GHOST: No ghost state for L = 3 | | | | siesta: System type = atom | Unit cell vectors (Ang): | 5.6699530.000.00 | 0.005.6699530.00 | 0.000.005.669953 | | | mulliken: Atomic and Orbital Populations: | | mulliken: Spin UP | | Species: C | Atom Qatom Qorb | 2s 2s 2py 2pz 2px 2py 2pz 2px | 2Pdxy 2Pdyz 2Pdz2 2Pdxz 2Pdx2-y2 | 1 3.000 1.013 -0.013 0.675 0.675 0.675 -0.008 -0.008 -0.008 | 0.000 0.000 0.000 0.000 0.000 | | mulliken: Qtot =3.000 | | mulliken: Spin DOWN | | Species: C | Atom Qatom Qorb | 2s 2s 2py 2pz 2px 2py 2pz 2px | 2Pdxy 2Pdyz 2Pdz2 2Pdxz 2Pdx2-y2 | 1 1.000 1.107 -0.107 0.000 0.000 0.000 0.000 0.000 0.000 | 0.000 0.000 0.000 0.000 0.000 | | mulliken: Qtot =1.000 | siesta: Total spin polarization (Qup-Qdown) =2.00 | | | siesta: Program's energy decomposition (eV): | siesta: Eions = 316.134481 | siesta: Ena = 109.046674 | siesta: Ekin=95.055041 | siesta: Enl = 6.977798 | siesta: DEna= 1.095781 | siesta: DUscf = 0.005992 | siesta: DUext = 0.00 | siesta: Exc = -43.532745 | siesta: eta*DQ = 0.00 | siesta: Emadel = 0.00 | siesta: Ekinion = 0.00 | siesta: Eharris = -147.485944 | siesta: Etot= -147.485941 | siesta: FreeEng = -147.535307 | | siesta: Final energy (eV): | siesta: Kinetic = 95.055041 | siesta: Hartree = 69.905598 | siesta:Ext. field = 0.00 | siesta: Exch.-corr. = -43.532745 | siesta: Ion-electron =-213.454182 | siesta: Ion-ion = -55.459652 | siesta: Ekinion = 0.00 | siesta: Total =-147.485941 | | siesta: Stress tensor (static) (eV/Ang**3): | siesta: 0.190.000.00 | siesta: 0.000.190.00 | siesta: 0.000.000.19 | | | | | | | | - | Do you Yahoo!? | Everyone is raving about the all-new Yahoo! Mail Beta.
Re: [SIESTA-L] -ve frequency in phonon calculation
On Mon, 5 Jun 2006, bipul rakshit wrote: | hi siesta users, | i am calculating the phonon frequencies of MgB2. i got nice match in bulk modulus of 156 GPa with LDA. While the experimental value is 151 GPa. | | With that Pseudopotential when i tried to calculate the phonon frequency by Vibra. i got -ve frequencies. | | can any body suggest, what can i do??? Bipul: You must obtain three acoustic modes with the frequencies close enough to zero, as the lowest eigenvalues. If they are ~ 0.01, either positive or negative, that's fine. If any of them largely deviates from zero, that is a indication that are not good because you have eggbox effect, as a result the acoustic sum rule is not satisfied with sufficient accuracy. You can check the eggbox effect prior to phonon calculation: displace all atoms uniformely, look at the variation of forces; the sum of forces over all atoms in each cartesian direction must be less than, at most, 0.1. Increase Mesh Cutoff to control the eggbox. Good equilibrium lattice constant and bulk modulus are good indication that the calculation is generally sane; however they are not so sensitive to the eggbox as the forces (and hence phonons) are. Good luck, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] how to calculate DOS of individual atom?
On Fri, 2 Jun 2006, Cherry Y. Yates wrote: | Dear all, | | I wonder if anyone know how to calculate electron DOS | of individual atom? I know there is a projected DOS block, but that one only records total DOS. For example, I have a system with 1000 atoms, I can calculate its total DOS, so that I can get band structure; but if I would like to calculate the electron DOS of the first 100 atoms, how could I do within SIESTA? Cherry: along with DOS you always get the .PDOS and pdos.xml files which contain electron DOS resolved in each single atom (and orbital). In order to select information in a sophisticated way (e.g., first 100 atoms as you wish), look in siesta-2.0/Util/pdosxml for inspiration. Good luck, Andrei Postnikov +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Phonon DOS
On Wed, 31 May 2006, Sharat Chandra wrote:
| Hi
| Can any one please tell me how to calculate the phonon DOS using Vibra
| module? Has any one developed a code for carrying out the calculation?
Hallo Sharat:
yes I have a tool to calculate (local) phonon DOS,
however it is of limited usefulness:
it assumes a single q-point in a presumable large supercell,
and smears out the peaks at frequencies (and with weights)
as provided by the vibrator.
Please find enclosed the source code and README.
Best regards,
Andrei
+-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+
| Paul Verlaine University - Institute de Physique Electronique et Chimie, |
| Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France |
+-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/
--+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 [EMAIL PROTECTED]
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 with a given Q-vector).'
101 write (6,703)
703 format(' Do you want total density of states (T) ',
. ' or convolution (Q) ? : ',$)
read (5,*,err=101,end=101) qlab
if (qlab.eq.'T'.or.qlab.eq.'t') then
mode=1
[Fwd: Re: [SIESTA-L] About dS12/dR12 calculated by matel]
On Thu, 26 Jan 2006, Daniel Sanchez Portal wrote: | Dear all, | | I have been checking a little bit the problem commented by | Andrei. From my point of view it seems a problem related to the use of a | interpolation from a (not very dense) grid. | We always have numerical values of S(R) and numerical derivatives. | The derivatives given by the splines seem indeed to be quite consistent | with the behaviour of given also for SR (I am talking now | about the radial part). | However, this behaviour | is not completely correct. If the number of points in the | radial mesh is increased one can see a significant | decrease in the magnitude of the spurious force. | | Daniel | Dear community, I agree with Daniel but I think this is a slightly different issue; Daniel writes about the numerical accuracy whereas my point was about missing symmetry. And this missing symmetry can be traced to the spherical harmonics part and has no relation whatsoever with the radial integration. I also checked that the error (deviation from zero of what must be zero) can be pushed down by increasing e.g. the k-cutoff in matel, but quite slowly. Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [Fwd: Re: [SIESTA-L] About dS12/dR12 calculated by matel]
Jose, Emilio: with respect to the matel issue discussed yesterday, I suggest now a more elegant substitution of my yesterday's fast and dirty fix. It explicitly uses the symmetry, that is that the spherical harmonics change with R - -R as (-1)**l whereas gradients as (-1)**(l-1). The changes in matel amount to two variables introduced and the final block re-written; they are all marked as C --- AVP begin . LVAL, C --- AVP end etc. in the attached matel.f I could reproduce my numerical tests of yesterday (and even better, with some noice of the order of ~E-29 now removed). However, my test system included only (s,p) functions so far. Please advise. I'll put this new fix in the mailing list as well, for the case someone would like to do more tests. Best regards, Andrei +-- Dr. Andrei Postnikov Tel. +33-387315873 - mobile +33-666784053 ---+ | Paul Verlaine University - Institute de Physique Electronique et Chimie, | | Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France | +-- [EMAIL PROTECTED] http://www.home.uni-osnabrueck.de/apostnik/ --+ ! ! Copyright (c) Fundacion General Universidad Autonoma de Madrid: ! E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal ! and J.M.Soler, 1996-2003 ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! SUBROUTINE MATEL( OPERAT, IS1, IS2, IO1, IO2, R12, S12, DSDR ) C *** C Finds the overlap or laplacian matrix elements between two atomic C orbitals. Written by J.M.Soler. April 1995. C Matrix elements of the position operator by DSP. June, 1999 C Last modifications by JMS. January 2001. C * INPUT *** C CHARACTER OPERAT : Operator to be used: 'S' = Unity (overlap) C 'T' = -Laplacian C 'X' = x C 'Y' = y C 'Z' = z C What is actually computed is the matrix elements: C 'X' =phi1(r-R12)|x|phi2(r) C 'Y' =phi1(r-R12)|y|phi2(r) C 'Z' =phi1(r-R12)|z|phi2(r) C i.e the origin is taken at the position of the second atom. C INTEGER IS1:Species index of 1st orbital atom C INTEGER IS2:Species index of 2nd orbital atom CAllowed range of IS1, IS2: (1, MAXS), where CMAXS and MAXL (below) are internal parameters. C INTEGER IO1: Orbital index of 1st orbital (within atom) C INTEGER IO2: Orbital index of 2nd orbital (within atom) CAllowed range of IO1, IO2: (-MAXL**2, MAXL**2) CIndexes IS1,IS2,IO1,IO2 are used to call routines CRCUT and PHIATM (see below). C REAL*8R12(3) : Vector from first to second atom C * OUTPUT ** C REAL*8 S12 : Matrix element between orbitals. C REAL*8 DSDR(3) : Derivative (gradient) of S12 with respect to R12. C * ROUTINES CALLED * C The following functions must exist: C C INTEGER FUNCTION LOFIO(IS,IO) C Returns the total angular momentum of orbitals and KB proyectors. C Input: C INTEGER IS : Species index C INTEGER IO : Orbital index C C REAL*8 FUNCTION RCUT(IS,IO) C Returns cutoff radius of orbitals and KB projectors. C Input: C INTEGER IS : Species index C INTEGER IO : Orbital index C C SUBROUTINE PHIATM(IS,IO,R,PHI,GRPHI) CFinds values and gradients of: Ca) basis orbitals (IO 0) Cb) KB proyectors (IO 0) Cc) Local part of screened pseudopotential (IO = 0) ( b) and c) are C not required if MATEL is called only with IO 0 ) C Input: C INTEGER IS : Species index C INTEGER IO : Orbital index C REAL*8 R(3) : Position with respect to atom C Output: C REAL*8 PHI : Value of orbital or KB projector at point R C REAL*8 GRPHI(3) : Gradient of PHI at point R C C SUBROUTINE XPHIATM(IS,IO,R,PHI,GRPHI) C The same as PHIATM but multiply by x. C C SUBROUTINE YPHIATM(IS,IO,R,PHI,GRPHI) C The same as PHIATM but multiply by y. C C
Re: [SIESTA-L] using Nose options
Dear Sergey, I think, the temperature in a simulation with Nose thermostat is supposed to fluctuate around the target value, because your system exchanges energy with the thermostat. Your fluctuataions seem quite wild, but probably that's because your test system is quite small. I guess they'll be damped in a larger system. Moreover there is the MD.NoseMass parameter responsible for coupling to the thermostat. However, I can give you no hints about how to tune its value. Consult the Nose paper (Progress of Theoretical Physics, Supplement No.103, 1991) - it is well readable. Best regards, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Fri, 15 Jul 2005, Sergey Lisenkov wrote: | Dear Siesters, | | I am a little confused when perform MD with Nose thermostat. I put TargetTemperature = 500 K but in output file I see the next behaviour: | | siesta: Temp_ion = 56.381 K | siesta: Temp_ion = 39.006 K | siesta: Temp_ion = 251.629 K | siesta: Temp_ion = 394.578 K | siesta: Temp_ion = 216.338 K | siesta: Temp_ion = 24.889 K | siesta: Temp_ion = 204.911 K | siesta: Temp_ion = 416.811 K | siesta: Temp_ion = 281.607 K | siesta: Temp_ion = 36.181 K | siesta: Temp_ion = 58.008 K | siesta: Temp_ion = 343.699 K | siesta: Temp_ion = 482.636 K | siesta: Temp_ion = 222.906 K | siesta: Temp_ion = 49.991 K | siesta: Temp_ion = 342.534 K | siesta: Temp_ion = 555.130 K | siesta: Temp_ion = 327.737 K | siesta: Temp_ion = 35.523 K | siesta: Temp_ion = 97.661 K | siesta: Temp_ion = 453.434 K | siesta: Temp_ion = 590.339 K | siesta: Temp_ion = 239.125 K | siesta: Temp_ion = 60.871 K | siesta: Temp_ion = 482.552 K | | So, a question is : how to perform MD calculations with constant temperature (Nose thermostat)? | | Many thanks, |Best wishes, | Sergey | |
[SIESTA-L] Conversion tools Siesta - XCrySDen
Dear Siesters, attached is a set of conversion tools, designed to transform Siesta properties files into XCrySDen input files (XSF for 2/3-dum. grids, BXSF for Fermi surface plotting, and AXSF for molecular dynamics vizualizations), see http://www.xcrysden.org/ Some of these tools were earlier subject to discussion on the List, but now I streamlined some things and made everything f90-conform. So whoever tried old ones (rho2xsf, eig2bxsf) please test new ones. I've done some preliminary tests on some systems, but sure enough there is always at least one more bug left. Enjoy and please send me back your criticisms and bug reports. Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ vers01.tar Description: Unix tar archive
[SIESTA-L] Fermi surface plotting
On Mon, 4 Jul 2005, xianghjun wrote: | Dear Andrei Postnikov, | I tried ene2bxsf, sometimes it works well, for example for factitial SCC Fe, | however ene2bxsf fails for BCC Fe, as an example in siesta. Dear H.J. Xiang, thank you for your hint. It was indeed a not clean enough piece of code which shoot for bcc lattice. Please try the updated version of eig2bxsf attached to this mail. Just one note of caution. In the Siesta example of Fe you use, the lattice vectors are 0.5 0.50 0.50 0.5 -0.50 0.50 0.5 0.50 -0.50 - this is of course legal for a band structure calculation, and even results in a Fermi surface, which is however weirdly folded and not quite of cubic symmetry. If you take more conventional and symmetric lattice vectors -0.5 0.50 0.50 0.5 -0.50 0.50 0.5 0.50 -0.50 the Fermi surface will be nice and cubic-symmetric. (However, I don't remember what the true Fe Fermi surface is like, so you better check it). Good luck, and again let me know in case of any problems. Andrei +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+C C ene2bxsf, a script to make the band-XSF file (for plotting C the Fermi surface) from siesta KP and EIG files. C The k-mesh must be generated without shift C (requirement of the BXSF format). C C Written by Andrei Postnikov, Jul 2005 Vers_0.1 C [EMAIL PROTECTED] C program ene2bxsf implicit none integer ii1,ii2,ii3,io1 parameter (ii1=11,ii2=12,ii3=13,io1=14) integer ii,jj,nbands,nbmin,nbmax,nkp,ikp,ndum,iis,mdiv(3),ind, .nspin,is,iband,iik,idiv1,idiv2,idiv3,itry1,itry2,itry3, .ib,mkp,homo(2),lumo(2),ik,in2,in3,ialloc character inpfil*60,outfil*60,syslab*30,suffix*6 real*8 cell(3,3),efermi,twopi,rcell(3,3),small,relmin(3), . kkp(3),dum real*8, allocatable :: relk(:,:), eneb(:,:), enek(:) integer, allocatable :: ndiv(:,:),irrek(:) parameter (small=1.0d-04) C C string manipulation functions in Fortran used below: C len_trim(string): returns the length of string C without trailing blank characters, C char(integer) : returns the character in the specified position C of computer's ASCII table, i.e. char(49)=1 C twopi = 8.d0*atan(1.d0) write (6,701) 701 format( Specify SystemLabel (or 'siesta' if none): ,$) read (5,*) syslab C --- open .XV and .EIG : inpfil = syslab(1:len_trim(syslab))//'.XV' open (ii1,file=inpfil,form='formatted',status='old',err=801) write (6,*) 'Found and opened: ',inpfil inpfil = syslab(1:len_trim(syslab))//'.EIG' open (ii2,file=inpfil,form='formatted',status='old',err=801) write (6,*) 'Found and opened: ',inpfil C --- read Fermi energy and total number of bands from .EIG : read (ii2,*) efermi read (ii2,*) nbands, nspin, nkp if (nspin.ne.1.and.nspin.ne.2) then write (6,*) 'A problem encountered: nspin=',nspin stop endif C --- finds bands crossing the Fermi energy: do is=1,nspin homo(is)=0 ! higest (partially?) occupied band lumo(is)=nbands+1! lowest (partially?) unoccupied band enddo write (6,*) ' nkp=',nkp,' nbands=',nbands allocate (eneb(nbands,nspin),STAT=ialloc) if (ialloc.ne.0) then write (6,*) ' Fails to allocate space for nbands=', .nbands,', nspin=',nspin stop endif do ik=1,nkp read (ii2,(i5,10f12.5,/,(5x,10f12.5))) iik, ! written in . ((eneb(ib,is),ib=1,nbands),is=1,min(nspin,2)) ! ioeif.f if (iik.ne.ik) then write (6,*) ' iik=',iik,'.ne. ik=',ik,' for spin ',is stop endif do is=1,nspin do ib=1,nbands if (eneb(ib,is).lt.efermi.and.ib.gt.homo(is)) homo(is) = ib if (eneb(ib,is).gt.efermi.and.ib.lt.lumo(is)) lumo(is) = ib enddo enddo enddo deallocate (eneb) do is=1,nspin write (6,*) ' is=',is,' homo, lumo=',homo(is),lumo(is) if (homo(is).lt.lumo(is)) write (6,201) is, homo(is),lumo(is) enddo C --- read cell vectors from .XV, convert to Ang, find reciprocal: do ii=1,3 read (ii1,*,end=803,err=803) (cell(ii,jj),jj=1,3) enddo close (ii1) call inver3(cell,rcell) C --- open .KP as old: inpfil = syslab(1:len_trim(syslab))//'.KP' open (ii3,file=inpfil,form='formatted',status='old',err=801) write (6,*) 'Found and opened: ',inpfil C --- read k-points from the .KP file and recover their fractional C coordinates in terms of reciprocal
Re: [SIESTA-L] grid2cube
On Sun, 26 Jun 2005, xianghjun wrote: | Dear Andrei Postnikov, | I have tried rho2xsf, and it works well. | Now I want to plot the Fermi surface using Xcrysden. | I learn that you have a program named ene2bxsf to do this job. | Would you please send a copy of ene2bxsf to me? | Thanks a lot. | | Best regards, | xianghjun Dear Xianghjun, please find in the attachment a code to calculate Fermi surface file in Xcrysden. I did not check it extensively, but it seems to work also with the version 1.4 (that with an additional option to plot several sheets of the Fermi surface in the same plot). Please let me know if you encounter any problems. Be sure that you use many enough k-points without shifting. Good luck, Andrei C C ene2bxsf, a script to make the band-XSF file (for plotting C the Fermi surface) from siesta KP and EIG files. C The k-mesh must be generated without shift C (requirement of the BXSF format). C C Written by Andrei Postnikov, Mar 2005 C [EMAIL PROTECTED] C program ene2bxsf implicit none integer MKP,MNB parameter (MKP=8000,MNB=500) integer ii1,ii2,ii3,io1 parameter (ii1=11,ii2=12,ii3=13,io1=14) integer ii,jj,nbands,nbmin,nbmax,nkp,ikp,ndum,iis, .mdiv1,mdiv2,mdiv3,ndiv(3,MKP),ind,irrek(MKP),nspin,is, .iband,iik,idiv1,idiv2,idiv3,itry1,itry2,itry3,ib, .homo(2),lumo(2),ik,in2,in3 character inpfil*60,outfil*60,syslab*30,suffix*6 real*8 cell(3,3),efermi,twopi,rcell(3,3),small,relmin(3), . relk(3,MKP),kkp(3),eneb(MNB,2),enek(MKP),dum equivalence (eneb,enek) parameter (small=1.0d-04) C C string manipulation functions in Fortran used below: C len_trim(string): returns the length of string C without trailing blank characters, C char(integer) : returns the character in the specified position C of computer's ASCII table, i.e. char(49)=1 C twopi = 8.d0*atan(1.d0) write (6,701) 701 format( Specify SystemLabel (or 'siesta' if none): ,$) read (5,*) syslab C --- open .XV and .EIG : inpfil = syslab(1:len_trim(syslab))//'.XV' open (ii1,file=inpfil,form='formatted',status='old',err=801) write (6,*) 'Found and opened: ',inpfil inpfil = syslab(1:len_trim(syslab))//'.EIG' open (ii2,file=inpfil,form='formatted',status='old',err=801) write (6,*) 'Found and opened: ',inpfil C --- read Fermi energy and total number of bands from .EIG : read (ii2,*) efermi read (ii2,*) nbands, nspin, nkp if (nbands.gt.MNB) then write (6,*) ' nbands=',nbands,' .gt. MNB=',MNB stop endif if (nspin.ne.1.and.nspin.ne.2) then write (6,*) 'A problem encountered: nspin=',nspin stop endif if (nkp.gt.MKP) then write (6,*) ' nkp=',nkp,' .gt. MKP=',MKP stop endif C --- finds bands crossing the Fermi energy: do is=1,nspin homo(is)=0 ! higest (partially?) occupied band lumo(is)=nbands+1! lowest (partially?) unoccupied band enddo write (6,*) ' nkp=',nkp,' nbands=',nbands do ik=1,nkp read (ii2,(i5,10f12.5,/,(5x,10f12.5))) iik, ! written in . ((eneb(ib,is),ib=1,nbands),is=1,min(nspin,2)) ! ioeif.f if (iik.ne.ik) then write (6,*) ' iik=',iik,'.ne. ik=',ik,' for spin ',is stop endif do is=1,nspin do ib=1,nbands if (eneb(ib,is).lt.efermi.and.ib.gt.homo(is)) homo(is) = ib if (eneb(ib,is).gt.efermi.and.ib.lt.lumo(is)) lumo(is) = ib enddo enddo enddo do is=1,nspin write (6,*) ' is=',is,' homo, lumo=',homo(is),lumo(is) if (homo(is).lt.lumo(is)) write (6,201) is, homo(is),lumo(is) enddo C --- read cell vectors from .XV, convert to Ang, find reciprocal: do ii=1,3 read (ii1,*,end=803,err=803) (cell(ii,jj),jj=1,3) enddo close (ii1) call inver3(cell,rcell) C write (6,*) ' cell:' C do ii=1,3 C write(6,'(3f15.6)') (cell(ii,jj),jj=1,3) C enddo C write (6,*) ' rcell:' C do ii=1,3 C write(6,'(3f15.6)') (rcell(ii,jj),jj=1,3) C enddo C --- open .KP as old: inpfil = syslab(1:len_trim(syslab))//'.KP' open (ii3,file=inpfil,form='formatted',status='old',err=801) write (6,*) 'Found and opened: ',inpfil C --- read k-points from the .KP file and recover their fractional C coordinates in terms of reciprocal vectors: read (ii3,*) nkp relmin(:)=1.d0/small do ikp=1,nkp read (ii3,*) iik,(kkp(jj),jj=1,3) if (iik.ne.ikp) then write (6,*) ' a mess in KP list: read in iik=',iik, .', but expected ikp=',ikp stop endif C --- find relative coordinates of k-point: do ii=1,3 relk(ii,ikp)=( cell(ii,1
Re: [SIESTA-L] diamond graphite phonon frequencies
Hi Mousumi, your zone-center acoustic frequency is not close enough to zero; this is an indication of bad, or badly converged, forces. (However you already get right sign and order of magnitude, that's not always automatically the case :-) Some considerations: 1. Your basis (SZ) is probably too small for a good phonon calculation, 2. your k-mesh (7*7*5) may be not fine enough. 3. Check if your MeshCutoff is sufficient to get small enough sum of forces over all atoms in X,Y,Z directions independently (200 Ry is not unreasonably low, but still...) Try first to get the bulk modulus in diamond right, this is not so lengthy calculation, so you can easily try different cutoffs. With a good bulk modulus, it's a fair chance that you'll get the phonons as well. Good luck, Andrei +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Sat, 21 May 2005, Mousumi Upadhyay Kahaly wrote: | Hi everybody, | | I am trying to relax then calculate the phonon frequencies of | different carbon systems using siesta fcbuild vibra. | | For diamond structure- at Gamma point: | my results- phonon freq= 0.3, 2042 cm-1. (only Gamma point calculation) | phonon freq= 0.50, 2163 cm-1 (Gamma to L-direction, 20 | points), | while expt result - 1333 cm-1. | | For graphite structure, at Gamma point, | my results- phonon freq.= 0.0, 1121.27 2305.7 cm-1, (only Gamma point | calculation) | | while expt result- 0.0, 900 1600 cm-1. | | Can anybody please suggest me, why are these discrepancies in | phonon frequencies? What should I do now to get rid of this | huge error? | | Here is my sample relaxation input file(that creats .FC file, | used to get phonon freq. through vibrator) for graphite -- | | best regards, mousumi | | *** | | SystemName Carbon in graphite # Descriptive name of the system | SystemLabel graphite # Short name for naming files | | NumberOfAtoms 4# Number of atoms | NumberOfSpecies 1# Number of species | | %block Chemical_Species_Label | 16C | %endblock Chemical_Species_Label | | PAO.BasisSize SZ | | # Lattice, coordinates, k-sampling | | %block kgrid_Monkhorst_Pack | 7 0 00.0 | 0 7 00.0 | 0 0 50.0 | %endblock kgrid_Monkhorst_Pack | | | AtomicCoordinatesFormat Fractional # Format for coordinates | AtomicCoorFormatOut Ang | | | %block AtomicCoordinatesAndAtomicSpecies |0.3334326744 0.3334326744 0.00 1 |0.6668653488 0.6668653488 0.00 1 |0.00 0.00 0.50 1 |0.3334326744 0.3334326744 0.50 1 | %endblock AtomicCoordinatesAndAtomicSpecies | | %block LatticeVectors | 1. 0. 0.00 | 0.5000 0.86602540378443864676 0.00 | 0.000.00 2.120 | %endblock LatticeVectors | | LatticeConstant2.510Ang | | | | #kgrid_cutoff7. Ang | | XC.functional LDA # Exchange-correlation functional type | XC.authors CA # Particular parametrization of xc func | SpinPolarized .false. # Spin unpolarized calculation | MeshCutoff 200. Ry # Equivalent planewave cutoff for the | grid | MaxSCFIterations100 # Maximum number of SCF iterations per | step | DM.MixingWeight 0.3 # New DM amount for next SCF cycle | DM.Tolerance1.d-4 # Tolerance in maximum difference | # between input and output DM | DM.NumberPulay 3 # Number of SCF steps between pulay | mixing | | # Eigenvalue problem: order-N or diagonalization | | SolutionMethod diagon # OrderN or Diagon | ElectronicTemperature 5 K# Temp. for Fermi smearing | | # Molecular dynamics and relaxations | | MD.TypeOfRunFC # Type of dynamics: | #ElectronicTemperature 30 K# Temp. for Fermi smearing | #MD.NumCGsteps 40 | | # Output options | | WriteCoorInitial | WriteCoorStep | WriteForces | WriteKpoints.true. | WriteEigenvalues.false. | WriteKbands .true. | WriteBands .true. | WriteMullikenPop1# Write Mulliken Population Analysis | WriteCoorXmol .false. | WriteMDCoorXmol .false. | WriteMDhistory .false. | WriteCoorXmol .false. | | # Options for saving/reading information | | DM.UseSaveDM
Re: [SIESTA-L] variation of lattice constants
On Tue, 10 May 2005, Prof.Dr. Kemal ÇOLAKOGLU wrote: | Dear A.Postnikov, | Many thanks for your suggestions on our problem.Firstly We have calculated | the lattice constant as 6.47 A (very cose to experiment , 6.48 A)for ground | state bulk properties.But we have set the input file more before as you | stated | MD.TargetPressure 0.5 GPa | MD.VariableCell T | to obtain the new equilibrium lattice constant for future work (Cij vs | P ).Our obtained lattice constants vs enegy | values are not expected ones as follows : | | lattice constantEnergy |6.0-1454,395571 |6.1-1454,390932 |6.2-1454,397768 |6.3-1454,397499 |6.4-1454,398580 |6.5-1454,395873 |6.6 - 1454,397538 |6.7 -1454,398341 |6.8-1454,382220 |6.9 -1454,385354 | As it can be seen , the data dont show the expected parabolic curve.So ,we | can not obtain new correct- lattice constant .What is the reason of this | discrepancy ? | Your hint will appreciated | Regards, | Kemal Colakoglu, | Gazi University, | Dept.of Phys. Dear Kemal, your data show indeed a tendency to fall around a parabola, however with huge fluctuations around, the reasons can be bad convergency? too few k-points? too low MeshCutoff? The usual suspects... In principle, you should check the stability of your results vs. k- and mesh cutoffs. Good luck, Andrei +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Problem in understanding .FC file
On Tue, 10 May 2005, Mousumi Upadhyay Kahaly wrote: | Dear All, | | I am trying to get phonon dispersion relation for carbon in | diamond structure. | | In order to get the dynamical matrix, I use MD.TypeOfRun FC | which prints out one .FC output file. | | The input structure has 2 atoms per unit cell and there are 3 | degrees of freedom. So I expect to get 6 * 6 dynamical matrix, ie | 36 elements. But in .FC file 72 elements are printed out. Could | anyone please help me to understand how the elements are written | in output .FC file and what do the 72 elements in my case mean? | Dear Mousumi, the .FC file contains not just the elements of dynamical matrix, but raw values of Force/Displacement, independently for each displacement which are all done in positive and negative direction. That's why there are 2* more elements than you expect. In the ideal case the corresponding lines must be identical, but in reality they are not, and the averaging over +/- helps to maintain accuracy. Specifically in your case: | Force constants matrix | 45.7520722 -0.0821671 -0.0826496 -- -F/d on atom #1 when displacing #1 along -X | -44.7156483 0.0799917 0.0801122 -- -F/d on atom #2 when displacing #1 along -X | 45.7417394 -0.0852655 -0.0851193 -- -F/d on atom #1 when displacing #1 along +X | -44.7073200 0.0833158 0.0830790 -- -F/d on atom #2 when displacing #1 along +X | -0.0822562 45.7518784 -0.0827168 -- -F/d on atom #1 when displacing #1 along -Y | 0.0798732-44.7154000 0.0802512 | -0.0849879 45.7414678 -0.0850387 | 0.0834938-44.7075596 0.0832066 | -0.0823887 -0.0825205 45.7516590 | 0.0798920 0.0800978-44.7152041 | -0.0851261 -0.0852046 45.7419141 | 0.0837935 0.0836693-44.7080396 | -44.7075637 0.0832532 0.0832499 -- -F/d on atom #1 when displacing #2 along -X | 45.0647213 -0.0855477 -0.0853743 -- -F/d on atom #2 when displacing #2 along -X | -44.7160376 0.0800102 0.0799784 | 45.0724592 -0.0837014 -0.0834930 | 0.0832718-44.7077354 0.0832417 | -0.0857780 45.0649619 -0.0854771 | 0.0800767-44.7161333 0.0799742 | -0.0833136 45.0722208 -0.0835637 | 0.0830177 0.0831134-44.7079753 | -0.0855199 -0.0854836 45.0650133 | 0.0798787 0.0799204-44.7159306 | -0.0831745 -0.0835189 45.0720666 -- -F/d on atom #2 when displacing #2 along +Z Mind that the .FC file contains -F/d in eV/(Ang^2). Best regards, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Treating of charged defects and pseudopotentials for ions
Many thanks to Riccardo Rurali and Emilio Artacho for their enlightening remarks about calculations with net charge! Now that this topic is under discussion, I wonder (if anybody has experience) whether we should check the convergency vs. cell size only when comparing DIFFERENTLY charged systems, or as well when we compare different cases regarding to the SAME net charge and number of atoms? Similarly: can we neglect the warning of not-SC, BCC, FCC and skip caring about the Madelung energy when comparing identically charged systems in identical cells? It seems that the Madelung correction from the background charge will be the same for these cases, independently on what we are doing with the molecule inside the cell, but the I may have overlooked something. Thank you ion advance, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Sat, 7 May 2005, Riccardo Rurali wrote: | Dear all, | | just a few minor remarks about the on-going discussion on charged systems. | | As far as I know, the neutralizing background charge is implemented, | regardless if the system is in cluster or supercell geometry. | As you know, this artificial background (needed so that the Coulomb energy | does not diverge) provokes a spourios interaction that is normally | (partially) corrected with the Madelung term. | | Now, what is true is that if your system is not a cluster the Madelung | correction is not calculated, even though the simulation box is of one of the | allowed type (SC, BCC or FCC). However, the reason is simply that, in the | case of a solid, the Madelung term needs to be scaled with the dielectric | constant and Siesta does not know it. | What I did in the past was commenting the appropiate IFs in the code in such a | way that the Madelung term is always calculated, but not summed to the total | energy. When the calculations was over, I read in output the value of the | Madelung term, scaled manually with the (experimental) dielectric constant | and add it to the total energy. | | Even so, if I got your point, you still have a problem: if your cell is not | SC, BCC or FCC you don't have the required alpha Madelung constant... | | A good idea could be calculating a large supercell. The obvious drawback with | respect to the option of a large cluster is that in this case the number of | atoms grows a lot (probably too much). In the limit of very large cell, the | amount of the Madelung correction (that you cannot calculate) becomes | negligible. | | Best Regards, | Riccardo | | -- | Riccardo Rurali | Laboratoire Collisions, Agrégats, Réactivité, IRSAMC | Université Paul Sabatier | 118 route de Narbonne | 31062 Toulouse | France | | e-mail: [EMAIL PROTECTED] | tel.+33 (0)5 61556071 | fax.+33 (0)5 61558317 | |
Re: [SIESTA-L] Treating of charged defects and pseudopotentials for ions
Hallo Konstantin, I think you mix up several things. The choice of pseudopot is one thing and a treatment of charged system is another. The pseudopot has nothing, or in a very indirect way if at all, to do with whether your system will be charged or not. In ideal case, the pseudopot is supposed to be transferable, i.e. able to treat cases with more or less electrons on a given atom on equal footing. If you are sure about a resulting charge state of your particular atom, you may wish to tune the pseudopot better for this case than the general one. But simply changing the electrion configuration in the pseudopot input does not make the resulting pseudopot suitable for a charge system. This will probably simply make it much much worse than before. To begin with, I won't care about tuning your pseudopot for Ni2+, I'd simply take a good one for Ni. Yours is not good. It is not a good idea to make spin-up and -down occupation numbers different in the pseudopot input unless you know what you are doing. That's because your resulting electron configuration of Ni impurity might be very far from 5(up)-3(down). Moreover your pseudoization radii are strange. (I did not test it though). My choice of Ni pseudopot is Ni pb rel pcec ATM3 no_date Troullier-Martins 4s 1.00r r= 2.00/4p 0.00r r= 2.00/3d 9.00r r= 2.00/4f 0.00r r= 1.50/ which probably can be improved in general, and/or tailored for your particular case. (If you care about magnetism, mind pseudocore correction). Then comes the charge of system. Whatever is provided in SIESTA is only for treating charged molecules or clusters, and don't be confused by the reference to SC, BCC and FCC in the corresponding section of the User's Guide - that's all about isolated impurities, put in simulation boxes of such shape. This is not your case. If you want to simulate charged impurities in a crystal, this won't help you because you won't have an extra charge in every cell of your crystal; the crystal will fall apart. What really happens in crystals with charged impurities is that your impurity gets or loses charge at the expense of some compensating center elsewhere in a crystal. The best you can do is to provide such charge-compensating mechanism in your calculation, or to perform a calculation with background charge (I am not sure this latter option is implemented). Finally you may be tempted to simulate a charged impurity in crystal by taking a (large) cluster as a model for your crystal, and charging the whole system, but I'm afraid this will create more problems than solve. Sorry for discouraging comments, but hopefully they might help anyway, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Fri, 6 May 2005, Konstantin Rushchanskii wrote: | Dear SIESTA users, | | I'm trying to investigate the influence of charged defects | on structural and electronic properties of some ferroelectric crystals. | For this purpose I need to construct a pseudopotential for the Ni +2 ion. | But I'm not sure in the way of choosing the valence and core states of | the ion. | So, if somebody could have a look at the input file, which I use, | and could check whether I'm doing the right thing or wrong, | it would be very appreciated. | | Here is my input file for neutral Ni atom: | | | pe Nickel, LDA, spin-polatized, rcore=1.06 |tm2 4.0 | n=Ni c=pbs | 0.0 0.0 0.0 0.0 0.0 0.0 |54 |40 1.00 1.00 |41 0.00 0.00 |32 5.00 3.00 |43 0.00 0.00 | 2 .80 3.15 0.79 2.90 1.00 1.06 | # | | And here is my input file for the Ni+2 ion: | | # | pe Nickel, LDA, spin-polatized, rcore=1.06 |tm2 4.0 | n=Ni c=pbs | 0.0 0.0 0.0 0.0 0.0 0.0 |54 |40 0.00 0.00 |41 0.00 0.00 |32 5.00 3.00 |43 0.00 0.00 | 2 .80 3.15 0.79 2.90 1.00 1.06 | ### | | I also would like to kindly ask, whether anybody in this community | has any experience in investigation of charged defects in non-cubic | supercells? | How should I take into account the neutralizing background in the case | of non-cubic supercell? | Is it real to obtain a converging result without using Madellung's | correction term, | which works only for cubic supercells? | | I will appreciate any help! | | Thank you very much in advance! | | -- | K.Rushchanskii, Dr. | Theoretische Physik, Universität Regensburg | 93040 Regensburg Germany | tel: +49 (0) 941
Re: [SIESTA-L] grid2cube
Dear Miguel and Rudy (and all),
now that we are discussing this topic, I also have a script
to convert RHO and other grid properties into XCrysden format.
I attach a Fortran source and a README file.
Your suggestions and criticisms are welcome.
(Miguel: you did not attach your script in your mail,
at least it didn't go to the list).
For whoever is interested, I also have a tool to plot Fermi surfaces
from SIESTA data with XCrysden.
As you may figure out, a tool to transfer wavefunctions with XCrysden is
not yet ready.
Best regards,
Andrei Postnikov
+-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+
| Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany |
+-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
On Mon, 25 Apr 2005, Miguel Pruneda wrote:
| Dear Rudy (and all),
|
| As far as I know, the format cube is only valid for orthogonal cells. There
| are other options to plot .RHO (and .DRHO, .VH, .VT, etc). Some of them are
| included in the siesta/Utils directory. Check the mailing-list history,
because
| this is a frequent topic of discussion.
|
| You can also use the format of XCrysden (called xfs) to plot 3D surfaces. I
| attach a script similar to grid2cube.f, written in f77, that you can use in
the
| same way as grid2cube. You can download XCrysden at: http://www.xcrysden.org/
|
| Cheers,
| Miguel
|
| Missatge citat per Rudy Coquet [EMAIL PROTECTED]:
|
| Dear SIESTA users,
|
| I would like to know if there is an easy way to convert a .RHO file to a
| cube file? Usually grid2cube works fine but I have a non-orthogonal
| cell so I cannot use it.
|
| Thank you,
|
| Rudy Coquet
| Cardiff Universtity, UK.
|
|
|
|
|
| -
| This mail sent through IMP: http://horde.org/imp/
|
|C
C rho2xsf, a script to transform 3-dim grid function
C (i.e. LDOS, RHO, DRHO, etc.) written in SIESTA
C by subr. iorho
C into a grid for XCrysden
C
C !!! IMPORTANT -- !!!
C compile this code with the same compiler switches as Siesta,
C otherwise reading the data from unformatted files will be spoiled.
C
C Written by Andrei Postnikov, Mar 2005
C [EMAIL PROTECTED]
C
program rho2xsf
implicit none
integer MPTS,mesh0(3),mesh1(3),ip,nspin,is,ii,jj,
.iat,nat,ityp,nz,ix,iy,iz,ind,mn,
.ixmax,iymax,izmax,ixmin,iymin,izmin
parameter (MPTS=1)
character inpfil*60,outfil*60,syslab*30,suffix*6
real*8 cell(3,3),coord(3),b2ang
parameter (b2ang=0.529177) ! Bohr to Angstroem
real func(MPTS),fdum ! NB! single precision, as in iorho.F
real fmax,fmin
C
C string manipulation functions in Fortran used below:
C len_trim(string): returns the length of string
C without trailing blank characters,
C char(integer) : returns the character in the specified position
C of computer's ASCII table, i.e. char(49)=1
C
write (6,701)
701 format( Specify SystemLabel (or 'siesta' if none): ,$)
read (5,*) syslab
inpfil = syslab(1:len_trim(syslab))//'.XV'
open (11,file=inpfil,form='formatted',status='old',err=801)
write (6,*) 'Found and opened: ',inpfil
C reads from .XV and writes into .XSF file in XCrysden format:
outfil = syslab(1:len_trim(syslab))//'.XSF'
open (12,file=outfil,form='formatted',status='new',err=802)
write (6,*) 'Opened as new:',outfil
C --- write crystal structure data:
write (12,'(a7)') 'CRYSTAL'
write (12,'(a7)') 'PRIMVEC'
do ii=1,3
read (11,*,end=803,err=803) (cell(ii,jj),jj=1,3)
write (12,204) (cell(ii,jj)*b2ang,jj=1,3)
enddo
write (12,'(a9)') 'PRIMCOORD'
read (11,*,end=804,err=804) nat
write (12,'(i4,a3)') nat,' 1'
do iat=1,nat
read (11,*,end=805,err=805) ityp, nz, (coord(ii),ii=1,3)
write (12,201) nz, (coord(ii)*b2ang,ii=1,3)
enddo
close (11)
C --- finished with .XV; now look for grid data files to include:
C
101 write (6,702)
702 format (' Add grid property (LDOS, RHO, ...;',
.' or BYE if none): ',$)
read (5,*) suffix
inpfil = syslab(1:len_trim(syslab))//
. '.'//suffix(1:len_trim(suffix))
open (11,file=inpfil,form='unformatted',status='old',err=806)
write (6,*) 'Found and opened: ',inpfil(1:len_trim(inpfil))
read (11,err=807) cell
read (11,err=808) mesh0, nspin
C
C --- Here, introdruce mesh step (keep every mn's point along
C each grid dimension).
11 write (6,'(a18,3i5)') 'Grid divisions: ',mesh0
write (6,'(a39,$)') Keep for XSF file every N'th point, N=
read (5,*) mn
if (mn.le.0) then
write (6,*) ' N must be positive, try again'
goto 11
endif
if (mod(mesh0(1),mn
Re: [SIESTA-L] SIESTA
On Tue, 5 Apr 2005, navaratnarajah kuganathan wrote: | I am Navaratnarajah Kuganathan am a DPhil Student working on Nanotubes | and Compounds using SIESTA in Oxford.I am now doing calculation on bulk | HgI2.I have used Fractional Coordinates according to the previous calculation | by researchers.I observed Density Matrix Converged to Zero for initial | (CG Move=0 ,within 50 SCF cycle).But it takes more than 200 for | CG move =1 ,2 ,3 etc. and difficult to converge. I am herewith attaching | the output file for my calculation. I would like to clarify | | 1.I have given correct fractional coordinates for HgI2. | For My electronic structure calculation, is it must to specify | CG steps = 50 or 100 etc.for geometry optimisation. | I think if we put correct values of Lattice parameters | and fractional coordinates,no need to optimise geometry. | | I would appreciate if you could help in this regard. Hallo, the answer depends on whether you really need to optimize geometry (then you must do it), or you want just electronic structure for a fixed lattice, correct values of Lattice parameters and fractional coordinates, in which case CG=0 is fine. Anyway, your calculation is probably wrong because you use AtomicCoordinatesFormat Fractional %block AtomicCoordinatesAndAtomicSpecies -1.0925 -1.0925 -3.11000 1 1.0925 1.09253.11000 1 ... that is weird (albeit legal) because fractional coordinates are supposed to be [ 0, 1 ]. Your cartesian coordinates (at the beginning of each CG move) are outcoor: Atomic coordinates (Ang): -4.77422500 -4.77422500 -38.68837388 1 Hg 1 4.774225004.77422500 38.68837388 1 Hg 2 ... that is probably not what you want, with your LatticeConstant4.370 Ang and LatticeParameters 1.001.002.84668 90. 90. 90. Moreover for structure relaxation in a crystal you'd need more k-points than one, and probably higher Mesh cutoff than default value. Good luck, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Generating WFS files
Hi, it is well hidden in the documentation but true: you should activate the block WaveFuncKPoints to make the program really write these coefficients. Good luck, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Fri, 4 Mar 2005, Fedwa El-Mellouhi wrote: | Hi, | | In the same context, I wanted Wavefuctions for the same tests. | Even if I specify WriteWaveFunctions .true., the program | doesn't provide any *.WFS files. | What's wrong ? | | Fedwa | | Le 4 mars 05, à 10:22, Mike F. Lisowski a écrit : | | Hello Everyone, | | I am trying to run Si in Siesta and want to get the wavefunction files | for | use with Denchar. I changed the longOutput .true. but I still do not | seem | to get the WFS files. What am I doing wrong? | | | *** | Fedwa El-Mellouhi tél: +1 (514) 343-6111 ext:1088 | Ph. D. Studentfax: +1 (514) 343-2071 | Département de physique | Université de Montréal | C.P. 6128, succ. Centre-ville | Montréal (Québec) H3C 3J7 | Canada[EMAIL PROTECTED] | *** | |
Re: [SIESTA-L] Mesh Supercell
On Wed, 9 Feb 2005, Aleksey Kletsov wrote: | Hello! | | How to cancel Internal auxiliary supercell calculation? I mean i want, | for example, to calculate Au2, so i do not need more atoms, just two. If you need just two (as an isolated cluster), you have to put them in large enough supercell, so that basis function won't overlap across the cell boundary. Good luck, Andrei Postnikov
Re: [SIESTA-L] Wave Functions: imaginary coefficients
On Tue, 1 Feb 2005, Aleksey Kletsov wrote: | Why some of wave functions do not have imaginary coefficients (they | are equal to zero)? Because those you are looking at, are at k=0. For other k values, they do. Good luck, Andrei Postnikov
Re: [SIESTA-L] Gold - Problem with factorization of B (part 2)
On Mon, 10 Jan 2005, Aleksey Kletsovfunky wrote: | Yes, they are default values. | | So what are the ways to avoid crashing diagonalization? To describe your structure correctly. At the moment you have two atoms in your cell which occupy identical positions. This is not what you want, and this is why your matrix becomes linear dependent, or whatever, and the diagonalization crashes. | P.S. What is this B? I don't know, nor have you to care about, once you fix your input. Good luck, Andrei +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] Gold - Problem with factorization of B
Hallo Aleksey, what are your lattice vectors? (check out.fdf) As you haven't define them, they are probably (100) (010) (001) by default. But you did define Lattice Constant of 2.5 Ang. The same is the distance between your two atoms in the cell. So they are sitting one on top of another, no wonder the diagonalization crashes. If you really need a dipole, put it in a large enough cell. Good luck, Andrei +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Mon, 10 Jan 2005, Aleksey Kletsovfunky wrote: | Thank you very much, Carlos! It worked! (i only changed number of | valence orbitals from 3 to 4) | But now i have another problem: | | The leading minor of order26 of B is not positive | definite. | The factorization of B could not be completed | and no eigenvalues or eigenvectors were computed. | Terminating due to failed diagonalisation | | Can anyone help me with this one? | | My .fdf file: | | SystemName Gold Dimer | SystemLabel Au2 | NumberOfAtoms 2 | NumberOfSpecies 1 | %block ChemicalSpeciesLabel | 1 79 Au | %endblock ChemicalSpeciesLabel | | %block PAO.Basis # Define Basis set |Au4 # Label, l_shells, type (opt), ionic_charge | (opt) | n=6 0 1 # n (opt if not using semicore levels), l, | Nzeta | 1.00 # rc(izeta=1,Nzeta)(Bohr) | n=6 1 2 | 1.00 1.00 | n=5 2 2 # l, Nzeta, PolOrb (opt), NzetaPol (opt) | 1.00 1.00 # rc(izeta=1,Nzeta)(Bohr) | n=5 3 1 | 1.00 | %endblock PAO.Basis | LatticeConstant 2.5 Ang | AtomicCoordinatesFormat Ang | %block AtomicCoordinatesAndAtomicSpecies | 0.000 0.000 0.00 1 | 0.000 0.000 2.50 1 | %endblock AtomicCoordinatesAndAtomicSpecies
Re: [SIESTA-L] One problem of using OrderN mothod.
Hi, it seems that OrderN works only for non-magnetic AND insulating systems, i.e. you must have i) an even total number of electrons (of course; as the simplest test that you are not trying to run a magnetic system), but also ii) your system must indeed be non-magnetic and iii) have a good band gap. Does your metal chain fulfil all these conditions? And: Do you really need OrderN for so small system? Best regards, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Tue, 4 Jan 2005, Xiao Changyong wrote: | Hi, all, | | When I chose the method OrderN to calculate a small system, single-cell | formed a metal wire, however, I got the following prompt: | | Ordern: enum= 99.000 | cspa: ERROR: Wrong total charge; odd charge: 0.E+000 | ERROR: Charge must be EVEN to use Order-N option. | | | In my calculation, a single-cell has 9 atoms, and each atom has 11 | valence electrons. | So that it means I must design a system with even number of electron. Am | I right? | If it is true, I must use a multi-cell for my case. | | I am looking forward to hearing from you. | | Regards, | David, Changyong XIAO |
Re: [SIESTA-L] Se problem
On Wed, 23 Jun 2004, Wenqing ZHANG wrote: | Hi, there, | | Does anyone have experience with Se bulk with Siesta? No, but with Se-based semiconductors. It seems essential to include 3d in the valence band for getting reasonable volume. I use Se ca nrl nc ATM3 no_date Troullier-Martins 4s 2.00 r= 1.89/4p 4.00 r= 1.99/3d10.00 r= 1.19/4f .00 r= 1.49/ and acordingly for the basis (DZP seems to work fine). Hopefully this helps... Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] comparison with abinit; core radius
On Mon, 14 Jun 2004, Jean-Paul Crocombette wrote: | Dear Raffaele, | | your Ni atoms are too close to get reasonnable results. | The pseudopotential spheres (defined by the maximum core radius) should | not overlap. | Your pseudopotential radius being 2.050 a.u. (1.1 Ang.), any calculation | with atoms closer than 2.2 Ang gives dubious results. | A distance of 1 Ang is much too small. Siesta probably complained about | the overlap in the output file. Dear Siesta developers: Is that true? Can anyone confirm or rebut this statement? While it must be so for a conventional pseudopot calculation, it seems I heard earlier that Siesta has a correction taking this into account, so that large core radii are no more a problem. With respect to Abinit, this might still be true (again, does anybody know for sure?). I think this is an important and sensitive issue. Thanks, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] How to optimize the bulk with SIESTA
On Thu, 13 May 2004, Hsin-Yen Chen wrote: | I like to optimize the BaTiO3 bulk in Cublic phase. | However, the MD.TypeOfRun CG was not converged with MD.NumCGsteps 500. | Furthermore, I checked that the output file to find the message, | siesta: System type = molecule. Did it be right? Not if you wanted to calculate bulk. Check your lattice parameters and internal coordinates. SIESTA thinks you have a molecule because your lattice parameters are too large. Cubic BaTiO3 has no internal degrees of freedom, so you would be better off just running a sequence of few volumes, in order to get optimized volume and bulk modulus, than doing a full CG search. | My object is to verify the BaTiO3 phonon properties and investigate the other | perovskite ferroelectics. Good luck. Check your electronic structure and bulk modulus against already available calculation data. Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] What is total energy?
On Tue, 2 Dec 2003, Sergey Lisenkov wrote: | DEar Prof. Soler, | | Thank you for your reply. Bu I cannot understand: Does SIESTA calculate the | real total energy? In most paper I see that authors compare total energy per | atom (for example, for carbon) with experimentasl results. How did they | obtain it? I checked it dor C-diamond, in the file carbpn.MDE the total | energy/per atom is more bigger than experimental results. | | Sergey | Dear Sergey, let's define what we speak about. The REAL total energy (i.e., reliable absolute values, like e.g. -1. Ry for hydrogen atom) is almost never available in DFT calculations, at least not with satisfactory accuracy. For one thing, the absolute numbers depend too much on implementation - limitations of the basis size, details of integrations etc. - and moreover there are basic limitations of DFT (spurious self interaction, that's why you'll never get -1 Ry for hydrogen; intraatomic correlations). One discusses absolute numbers in many-body calculations like configuration interaction, but then, for really pushing to the limits, you are technically restricted to one or few atoms. The absolute total energy in Siesta are wrong even by the order of magnitude, because the code uses pseudopotentials, so the (large) core contribution to the total energy is not present. However, this is not so bad, because if you are after ABSOLUTE total energy of the carbon atom SIESTA won't be the code of your choice. The RELATIVE properties like cohesion energies, defect formation energies etc. usually come about reasonably in density functional calculations, and also (with a bit of luck) with SIESTA. That's because these values are difference between two total energies, and the major part of systematic error cancels (including the missing core contribution). Nevertheless, you have to be careful to check that these difference converge depending on your basis size, energy cutoff etc. Hopefully this helps, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+
Re: [SIESTA-L] How to read the *.PDOS file
In the attachment there is a short Fortran program I wrote to manipulate PDOS files. Good luck, Andrei Postnikov +-- Dr. habil. Andrei Postnikov - Tel. +49-541-969.2377 -- Fax .2351 ---+ | Universitaet Osnabrueck - Fachbereich Physik, D-49069 Osnabrueck, Germany | +-- [EMAIL PROTECTED] - http://www.home.uni-osnabrueck.de/apostnik/ --+ On Mon, 17 Nov 2003 [EMAIL PROTECTED] wrote: | Dear siesta users and developers: | Who can give me some advice on how to read the *.PDOS file. I want to plot the | projected DOS,however, I don't know What kind of plot tool can use to read it and | then plot. Should I write a program to read it and then plot by myself? | Thanks and Best wishes! | Attachment Converted: c:\cygwin\home\paz\e-mail\attach\fmpdos.f fmpdos.f Description: Binary data
