[Wien] mbj on Diamond

2010-10-15 Thread Peter Blaha
Hi,
After I got the struct file, I could verify the problem.

As expected, it is again in the SRC_lapw0/brj.f subroutine,
where one has to find the root of a function.

For strange densities this is numerically non-trivial.

The problem at the nucleus of heavy elements was solved before, but
here is the problem in the interstitial, when rho is very small and
also grad-rho, tau and laplace-rho are sufficiently small.

Then the required functions are nearly zero (lt. 10^-6) for a range
of x=10-30; but using x=30 produces a V-xc potential of -100 Ry,
which is the reason for your Eigenvalues below zero.

When such problems occur again, please check also case.output0.
The Fourier coefficients of Vxc must converge, i.e. (0 0 0) should be
order one, while (0 0 30) should be order 10^-5 .

The attached subroutine brj.f should fix these problems (at least your case 
converges
smoothly).

Dear All,



We are performing the mbj calculations for a carbon based compound. According 
to the usersguide there are three SCF cycles for mbj calculations: first a 
regular calculations 
within LDA/GGA (we use the PBE-GGA here),  second one more cycle run_lapw ?NI 
?i 1 , and third the mbj run after changing the potential energy functional 
indxc=28 in case.in0 and 
index=50 in case.in0_grr.

Here we call the regular SCF cycles C1.scf, second one-more SCF cycle as 
C2.scf, and the third the mbj as cycle C3.scf.

The first regular cycle and the second run_lapw ?NI ?i 1 are converged 
smoothly. However, the third mbj cycle is stopped at lapw2 in its second 
iteration.

We analyzed the problem to find the source of the error. The result is given 
below, where the C2.scf line refers to the last :ITE of the second one more SCF 
cycle, and the C3.scf 
refers to the first :ITE of the third mbj run:



C2.scf::NTO033: TOTAL   CHARGE IN SPHERE  1 =3.9781366

C3.scf::NTO033: TOTAL   CHARGE IN SPHERE  1 =2.4250427



C2.scf::CTO033: TOTAL   CHARGE IN SPHERE  1 =3.9781254

C3.scf::CTO033: TOTAL   CHARGE IN SPHERE  1 =3.9470631



C2.scf::DIS  :  CHARGE DISTANCE   ( 0.355 for atom   33 spin 1)  
0.136

C3.scf::DIS  :  CHARGE DISTANCE   ( 1.8978668 for atom   25 spin 1)  
1.5016586



C2.scf::NEC01: NUCLEAR AND ELECTRONIC CHARGE366.0   365.98257 
1.5

C3.scf::NEC01: NUCLEAR AND ELECTRONIC CHARGE366.0   365.98171 
1.5



C2.scf::FER  : F E R M I - ENERGY(TETRAH.M.)=   0.21390

C3.scf::FER  : F E R M I - ENERGY(TETRAH.M.)=  -1.44751



The result clearly shows that there is a jump in :NTO, :DIS, and :FER (But in 
:CTO) after changing the functional to index=28.



-- 

   P.Blaha
--
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-15671 FAX: +43-1-58801-15698
Email: blaha at theochem.tuwien.ac.atWWW: http://info.tuwien.ac.at/theochem/
--
-- next part --
A non-text attachment was scrubbed...
Name: brj.f
Type: text/x-fortran
Size: 3680 bytes
Desc: not available
URL: 
http://zeus.theochem.tuwien.ac.at/pipermail/wien/attachments/20101015/f8063dbe/attachment.f


[Wien] MAE of Co hcp

2010-10-15 Thread Peter Blaha
No, you are NOT using P1 symmetry (but you have 24 symops.)

As I said, you have to find a common symmetry for all your cases. Your struct 
file definitely
does not do this.

If you are unable to manage this by yourself, eventually you may have to use a 
trick like
adding a 3rd atom at some arbitrary position, so that the WIEN2k 
initialization programs
really find P1 symmetry with only ONE symop. (or the common one). Stop the 
initialization
after kgen.
After you have this, remove the 3rd atom from struct and in* files and run 
lstart and dstart
to complete the initialization.
The runsp
save_lapw
The do the spinorbit+lapw2 for the two directions.

PS: Please note: As far as I remember, the MAE of Co is only a few mycroRy, 
i.e. one had to use
enormous k-meshes to get anything meaningful (and LDA/GGA is wrong anyway ?). 
Check literature.


Am 13.10.2010 10:27, schrieb Bin Shao:
 Dear Prof. Peter Blaha,

 With your suggestion, I recalculate the hcp Co in P1 symmetry with experiment 
 parameters, and the k-mesh is 40x40x24. Then I get the energy in (001), (100) 
 and (010) direction.
 But I find that the energy in (100) is large than (001) and that in (010) is 
 smaller than (001). In my opinion, they should be the same in (100) and in 
 (010) direction.

 Please give me some comments, and the attachment is my struct file.

 Thank you in advance!

 Best,

 On Tue, Oct 12, 2010 at 2:46 PM, Peter Blaha pblaha at theochem.tuwien.ac.at 
 mailto:pblaha at theochem.tuwien.ac.at wrote:

 If you change the k-mesh, you need to recalculate also lapw1.

 However, the proper approach is to find a symmetry (usually the lower 
 one), which accommodates
 both directions of the magnetization (usually you can run also the 
 higher symmetry case with the
 low symmetry). Run already the non-SO calculation in this symmetry and 
 then use only
 lapwso/lapw2 with the two directions in case.inso.



 Am 12.10.2010 05:24, schrieb Bin Shao:

 Dear all,

 I intend to calculate the MAE of hcp Co with force theorem. After the 
 nosoc scf calculation, I add the soc non-selfconsistently with the directions 
 of M || c and M || a,
 respectively. I use the initso_lapw to prepare the input files which 
 creats new structs and new klists for spin-polarized case. Then I run the 
 program by commands

 x lapwso -up -p
 x lapw2 -c -up -p
 x lapw2 -c -dn -p

 But here comes some errors in lapw2,

 dnlapw2.error
 
 
 'FERMI' - number of k-points inconsistent when reading kgen
 'FERMI' - check IN1 and KGEN files!
 **  testerror: Error in Parallel LAPW2
 
 

 Before Co hcp, I calculated the case of Fe monolayer with the same 
 approach, but no errors. How to deal with the problem or do I need to provide 
 some more input files?

 Any suggestion will be appriciated, thank you in advance!

 Best regards,

 --
 Bin Shao, Ph.D. Candidate
 College of Information Technical Science, Nankai University
 94 Weijin Rd. Nankai Dist. Tianjin 300071, China
 Email: binshao1118 at gmail.com mailto:binshao1118 at gmail.com 
 mailto:binshao1118 at gmail.com mailto:binshao1118 at gmail.com



 ___
 Wien mailing list
 Wien at zeus.theochem.tuwien.ac.at mailto:Wien at 
 zeus.theochem.tuwien.ac.at
 http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien


 --

   P.Blaha
 --
 Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
 Phone: +43-1-58801-15671 FAX: +43-1-58801-15698
 Email: blaha at theochem.tuwien.ac.at mailto:blaha at 
 theochem.tuwien.ac.atWWW: http://info.tuwien.ac.at/theochem/
 --
 ___
 Wien mailing list
 Wien at zeus.theochem.tuwien.ac.at mailto:Wien at 
 zeus.theochem.tuwien.ac.at
 http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien




 --
 Bin Shao, Ph.D. Candidate
 College of Information Technical Science, Nankai University
 94 Weijin Rd. Nankai Dist. Tianjin 300071, China
 Email: binshao1118 at gmail.com mailto:binshao1118 at gmail.com



 ___
 Wien mailing list
 Wien at zeus.theochem.tuwien.ac.at
 http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien

-- 

   P.Blaha
--
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-15671 

[Wien] optical properties calculations?

2010-10-15 Thread Jian-Xin Zhu
Dear Peter,

Thank you very much for the clarification. I have some comments below
and your further inputs are greatly appreciated.

On Oct 14, 2010, at 11:44 PM, Peter Blaha wrote:

 In the UG, page 147, I notice the following for the optical properties
 calculations 
 
 In cases of non-spinpolarized calculations WITHOUT inversion symmetry
 AND spin-orbit coupling,
 one must do some tricks and ?mimick? a spinpolarized calculation:
 
 Sorry for the bad English. It should read:
 In cases of non-spinpolarized spin-orbit calculations WITHOUT inversion 
 symmetry 
 
 As I understand, the first sentence is saying that the case is non-spin
 polarized but has no Inversion symmetry and no spin-orbit coupling.
 

For the cases of non-spinpolarized spin-orbit coupling, I can also take the 
following procedure to do the calculations ---

First change TOT to FERMI  and also use TETRA with a value of 101.0 in case.in2c

   (run_lapw) options: -so -s lapw1 -e lcore -p
Thu Oct 14 21:44:00 MDT 2010 (x) lapw1 -p
Thu Oct 14 21:46:18 MDT 2010 (x) lapwso -p
Thu Oct 14 21:49:23 MDT 2010 (x) lapw2 -c -so -p
Thu Oct 14 21:49:28 MDT 2010 (x) lcore
Thu Oct 14 21:58:23 MDT 2010 (x) opticc -so -p
Thu Oct 14 22:05:58 MDT 2010 (x) joint
Thu Oct 14 22:30:37 MDT 2010 (x) kram

Is there anything wrong with my procedure?

Can I understand the purpose of the tricks mentioned in UG is to   mimick a 
spin polarized calculation ?

BTW, I did not use p-1/2 relaticvistic LOs in LAPWSO as warned in the UG. 

Does the current version of OPTICS now support  p-1/2 relativistic LOs?


 2/ Another question:
 
 Does the current OPTICS support the LDA+U?
 
 Yes.

That's great. Will the following procedure do the job when a spin-polarized 
spin-coupling LDA+U case is considered? 

 
First change TOT to FERMI  and also use TETRA with a value of 101.0 in case.in2c

runsp_lapw -so -orb -s lapw1 -e lcore
x opticc -so -up
x joint -up
x kram -up 

Please note I only add the -orb option in the line 
runsp_lapw -so -orb -s lapw1 -e lcore

not in the other lines. 

For the forced non-spin polarized spin-orbit coupling LDA+U case, I 
would simply replace the runsp_lapw by runsp_c_lapw. Does it make sense?


Thank you very much for the instruction. 

Jianxin



 -- 
 Peter Blaha
 Inst.Materialchemie, TU Wien
 Getreidemarkt 9
 A-1060 Vienna
 Austria
 ___
 Wien mailing list
 Wien at zeus.theochem.tuwien.ac.at
 http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien

-- 
###
Jian-Xin Zhu, Ph.D
Theoretical Division, MS B262
Los Alamos National Laboratory
Los Alamos, New Mexico 87545
Phone: (505) 667 2363
Fax: (505) 665 4063
Email (main): jxzhu at lanl.gov
Email (backup): physjxzhu at gmail.com
URL: http://theory.lanl.gov
###








[Wien] Wien2k: Rhombohedral lattice

2010-10-15 Thread Jian-Xin Zhu
Dear Wien Users, 

This is an earlier communication I made with Peter.
It might be useful to be documented in the archive.

Regards, 

Jianxin

Begin forwarded message:

 From: Peter Blaha pblaha at theochem.tuwien.ac.at
 Date: September 18, 2010 1:31:45 AM MDT
 To: Jian-Xin Zhu jxzhu at lanl.gov
 Subject: Re: Wien2k: Rhombohedral lattice
 
 Yes, everithing is as statet in the UG (even if it is very unlogical
 to specify hexagonal lattice parameters and rhombohedral coordinates).
 
 PS: This is an IDEAL question for the mailing list. Others could learn
 from the answer too and I do not have to  answer this 15 times/year.
 
 Jian-Xin Zhu schrieb:
 Dear Peter, I am trying to calculate a system of rhombohedral lattice.  From 
 the UG, I  am advised that special care should be taken for preparing the 
 structure file, e.g., rhomb.struct. ---
 1/ On page 39,  for the input of unit cell parameters (a,b,c), it reads ...
 for rhombohedral (R) lattices the hexagonal lattice constants must be 
 specified. (The following may help you to convert between hexagonal and 
 rhombohedral specifications:
  $a_{hex} = 2 cos (\frac{\pi- \alpha_{rhomb}}{2} ) a_{rhomb}$
  $c_{hex} = 3 \sqrt{a_{rhomb}^2 - \frac{1}{3} a_{hex}^2 } $
  and (for fcc-like lattices) $a_{rhomb}=a_{cubic}/\sqrt{2} $
 
 It means I should put in the lattice parameter of the corresponding 
 hexagonal lattice. So in the above formula, alpha_rhomb is the angle between 
 any of two rhombohedral  primitive unit vectors, for example, alpha_homb=60 
 degree. Is it correct?
 2/ Then for the position of atom in internal units (x,y,z), it reads ...
 For R lattice use rhombohedral coordinates. (To convert from hexagonal into 
 rhombohedral coordinates use the auxiliary program *hex2rhomb*, which can be 
 called at a command-line:
  $ \vec X_{ortho} = \vec X_{hex} \left ( \begin{array}{ccc} 0  
 1  0 \\ \frac{\sqrt{3}}{2}  \frac{-1}{2}  0 \\ 0  0  1 \end{array} 
 \right ) $
  $ \vec X_{rhomb} = \vec X_{ortho} \left ( \begin{array}{ccc} 
 \frac{1}{\sqrt{3}}... ...rt{3}}  \frac{-2}{\sqrt{3}}\\ -1 1  0 \\ 1  1  
 1 \end{array} \right ) $
 
 Here when you say For R lattice use rhombohedral coordinates. Does it 
 means I should use the following x,y,z values  in
 vecr = x veca_rhombohedral + y vecb_rhombohedral + z vecc_rhombohedal
 rather than in
 vecr = x veca_hexagonal + y vecb_hexagonal + z vecc_hexagonal
 although we use the hexagonal lattice constants for the unit cell parameters 
 of the rhombohedral lattice? Here veca_rhombohedral, vecb_rhombohedral, 
 vecc_rhombohedral, are three primitive translation vectors for the 
 rhombohedral lattice I am going to study while veca_hexagonal, 
 vecb_hexagonal, vecc_hexagonal, are those for the hexagonal lattice 
 mentioned in the above point 1.
 I want to make sure I really understand these points of care correctly.
 Best regards, Jian-Xin
 P.S.: Since it is not about the code itself and also might be too trivial 
 for other users, I don't post it to the mailing list.   --
 ###
 Jian-Xin Zhu, Ph.D
 Theorertical Division, MS B262
 Los Alamos National Laboratory
 Los Alamos, NM 87545
 Phone: (505) 667 2363
 Fax: (505) 665 4063
 Emai: jxzhu at lanl.gov mailto:jxzhu at lanl.gov
 Email (backup): physjxzhu at gmail.com mailto:physjxzhu at gmail.com
 URL: http://theory.lanl.gov
 ###
 
 -- 
 -
 Peter Blaha
 Inst. Materials Chemistry, TU Vienna
 Getreidemarkt 9, A-1060 Vienna, Austria
 Tel: +43-1-5880115671
 Fax: +43-1-5880115698
 email: pblaha at theochem.tuwien.ac.at
 -

-- 
###
Jian-Xin Zhu, Ph.D
Theoretical Division, MS B262
Los Alamos National Laboratory
Los Alamos, New Mexico 87545
Phone: (505) 667 2363
Fax: (505) 665 4063
Email (main): jxzhu at lanl.gov
Email (backup): physjxzhu at gmail.com
URL: http://theory.lanl.gov
###






-- next part --
An HTML attachment was scrubbed...
URL: 
http://zeus.theochem.tuwien.ac.at/pipermail/wien/attachments/20101015/4a80bcb8/attachment.htm