Hi Andrie, A couple of comments, first the absolute value of the polarization per se is only defined modulo a quantum of polarization, so the value you get is P = P_"real" + n P_q, where n is an integer number and P_q is the quantum of polarization. P_q=e/A where A is the surface area of your unit cell perpendicular to the direction with respect to which you are calculating P.
Since you don't know the value of n, to get the value of P_"real" you usually calculate the polarization of a reference structure which P is known by symmetry (in perovskites for instance, the centrosymetric phase has either 0 or P_q/2). Then substracting the P from the centrosymmetric you should be able to get the polarization of the polar phase (you might still have some problems to identify the "true" polarization if the P of the polar phase is comparable or larger than the quantum, see this introductory paper for more details http://arxiv.org/abs/1202.1831) Another issue is the PolarizationGrid block you are using, which is the example in the siesta manual and might not be the most suitable for you simulation cell. Something like this would make more sense (re-read the manual entry for the meaning of this block): %block PolarizationGrids 20 4 4 yes 4 20 4 yes 4 4 15 yes %endblock PolarizationGrids Best regards, Pablo On Tue, Oct 1, 2013 at 2:39 PM, Andrei Buin <[email protected]> wrote: > Dear Siesta forum, > > I have 48(pervoskite CH3NH3PbI3) atoms in the 9 Angs x 9 Angs x 12 Angs > unit cell. > I'm trying to compute the Berry phase polarization usign > PolarizationGrids, and i get > insane dipole moment of: > > siesta: Macroscopic polarization per unit cell (Debye): > siesta: Along the lattice vectors 2742.652644 2679.284584 4191.471845 > siesta: Along cartesian directions 2742.652644 2679.284584 4191.471845 > > > > Input is attached. Forces are already converged. > > # 6.2 General System descriptors > > SystemName Tetra Super Unit Cell Siesta Calculations > SystemLabel Tet # Short name for naming files > NumberOfSpecies 5 > NumberOfAtoms 48 > > %block Chemical_Species_Label > 1 82 Pb > 2 53 I > 3 7 N > 4 6 C > 5 1 H > %endblock Chemical_Species_Label > > > > > # 6.3 Basis definitions > > PAO.BasisSize DZP > PAO.EnergyShift 65 meV > #PAO.SplitNorm 0.15 > PAO.SplitTailNorm true > PAO.SoftDefault true > PAO.OldStylePolOrbs false > > %block PS.KBprojectors > I 4 > 0 2 > 2000 -2000 > 1 2 > 2000 -2000 > 2 2 > 2000 -2000 > 3 1 > 2000 > > > Pb 4 > 0 2 > 2000 -2000 > 1 2 > 2000 -2000 > 2 2 > 2000 -2000 > 3 1 > 2000 > > N 4 > 0 2 > 2000 -2000 > 1 2 > 2000 -2000 > 2 2 > 2000 -2000 > 3 1 > 2000 > > H 3 > 0 2 > 2000 -2000 > 1 2 > 2000 -2000 > 2 1 > 2000 > > C 4 > 0 2 > 2000 -2000 > 1 2 > 2000 -2000 > 2 2 > 2000 -2000 > 3 1 > 2000 > %endblock PS.KBprojectors > > > > %block PAO.Basis # Define Basis set > Pb 3 # Species label, number of l-shells > n=6 0 2 # n, l, Nzeta, Polarization, NzetaPol > 0 0 > n=6 1 2 P 1 # n, l, Nzeta, Polarization, > NzetaPol > 0 0 > n=5 2 1 # SZ for 5d > 0 #0 > > > # H in C6H6 DZP > #Vova > H 2 0.00000 > n=1 0 2 E 11.36136 0.00928 > 7.72405 2.19949 > n=2 1 1 E 41.15301 0.00947 > 2.89938 > > > # C in C6H6 DZP > #Vova > C 3 0.00000 > n=2 0 2 E 39.65304 6.21693 > 7.40483 4.90026 > n=2 1 2 E 27.05294 3.74121 > 7.88345 3.11808 > n=3 2 1 E 55.60264 0.01540 > 3.93573 > %endblock PAO.Basis > > > > > > # 6.4 Lattice, coordinates, k-sampling > > LatticeConstant 8.96 Ang # 6.05 350Ry, 5.93845 Ang Exp, 5.936 Exp, > 5.84_LDA_Zhenya > > %block LatticeVectors > 1.0 0 0 #18A crystal + 12Avacuum > 0 1.0 0 > 0 0 1.4375 > %endblock LatticeVectors > > AtomicCoordinatesFormat Ang > AtomicCoorFormatOut Ang > > > > > > > > > > %block AtomicCoordinatesAndAtomicSpecies > -0.03167 -0.11487 -0.35971 1 > -0.12412 -0.42211 2.92565 2 > 1.85899 2.43724 -0.06194 2 > 4.40361 4.31212 -0.38969 1 > 6.99200 6.23198 -0.88879 2 > 2.46456 6.94775 -0.15215 2 > 6.36927 1.84289 0.08180 2 > -0.09247 -0.14207 6.12437 1 > 2.57611 1.64461 5.69084 2 > 4.35503 4.31005 6.09334 1 > 6.18152 6.97008 6.46046 2 > 4.59235 4.49865 2.91059 2 > 4.43687 3.75124 9.32494 2 > 0.44143 -0.12135 9.34642 2 > 3.78065 0.42371 2.39084 3 > 4.82458 -0.31863 3.15069 4 > 3.28270 -0.20569 1.72441 5 > 4.20292 1.18515 1.82761 5 > 3.08493 0.84035 3.03696 5 > 5.25142 0.34763 3.90531 5 > 5.60253 -0.63714 2.45202 5 > 4.36520 -1.18663 3.63088 5 > 4.08008 0.15264 8.72657 3 > 4.82907 -0.48545 9.84708 4 > 4.39959 -0.21875 7.81266 5 > 4.21776 1.18640 8.72722 5 > 3.05391 -0.01415 8.80689 5 > 4.48385 -0.06168 10.79380 5 > 5.89626 -0.28577 9.71242 5 > 4.64461 -1.56270 9.82703 5 > -0.85527 5.05646 2.32772 3 > 0.08950 4.13502 3.02005 4 > -0.81071 4.97580 1.29648 5 > -0.65212 6.05774 2.55366 5 > -1.84958 4.87012 2.59300 5 > -0.01489 4.27318 4.09978 5 > 1.10925 4.37398 2.70535 5 > -0.15654 3.10553 2.74484 5 > 7.04591 2.46318 6.25281 2 > 1.72806 6.16802 6.24656 2 > -0.59017 4.78536 8.67712 3 > 0.30844 3.96725 9.53987 4 > -1.16287 4.17787 8.05170 5 > -0.03405 5.40998 8.05433 5 > -1.22474 5.36318 9.25668 5 > 0.88403 4.63525 10.18598 5 > 0.97946 3.39039 8.89784 5 > -0.30006 3.29480 10.15037 5 > %endblock AtomicCoordinatesAndAtomicSpecies > > > #%block GeometryConstraints > #position 1 2 > #%endblock GeometryConstraints > > #kgrid_cutoff 20 Ang > > #BandLinesScale pi/a > #BandLinesScale ReciprocalLatticeVectors > > #%block BandLines > #1 1 1 1 L > #20 0 0 0 G > #20 1.5 0 1.5 K > #10 2 0 0 X > #20 0 0 0 G > #%endblock Bandlines > > %block kgrid_Monkhorst_Pack > 8 0 0 0 > 0 8 0 0 > 0 0 8 0 > %endblock kgrid_Monkhorst_Pack > > > > > # 6.5 DFT, Grid, SCF > > XC.functional GGA > XC.authors PBE > #SpinPolarized true > MeshCutoff 375 Ry > FilterCutoff 375 Ry > MaxSCFIterations 199 # Maximum number of SCF iter > DM.MixingWeight 0.05 # New DM amount for next SCF cycle > DM.NumberPulay 10 > DM.PulayOnFile false > DM.NumberKick 50 > #DM.KickMixingWeight 0.1 > #DM.Tolerance 0.00001 # 0.0004 > #DM.EnergyTolerance 0.00005 eV > > > > > > #%block PolarizationGrids > # 7 3 3 yes > # 3 7 3 yes > # 3 3 7 yes > #%endblock PolarizationGrids > > > %block PolarizationGrids > 10 3 4 yes > 2 20 2 no > 4 4 15 > %endblock PolarizationGrids > > > > # 6.6 Eigenvalue problem: order-N or diagonalization > > SolutionMethod diagon > #Diag.DivideAndConquer false > #NumberOfEigenStates 10000 # total is ~27000 SZ > #OccupationFunction MP > #OccupationMPOrder 1 > #ElectronicTemperature 300 K > #ON.ChemicalPotentialUse true > > > > # 6.7 Molecular dynamics and relaxations > > MD.TypeOfRun CG > #MD.NoseMass 600 Ry*fs**2 > > MD.NumCGsteps 9950 > MD.MaxCGDispl 0.15 Ang > MD.MaxForceTol 0.040 eV/Ang > > #MD.LengthTimeStep 1 fs > #MD.FinalTimeStep 10000 > #MD.InitialTemperature 550 K > #MD.TargetTemperature 550 K > > > > > > > # 6.8 Parallel options > > #BlockSize 426 # for 16 procs #852 for 8 # DZP 6813 orbitals > #ProcessorY 8 > #Diag.Memory 4 > #Diag.ParallelOverK true > > > > # 6.9 Efficiency options > > # 6.10 Output options > > WriteMDXmol true > WriteCoorXmol true > WriteMDhistory true > WriteXML false > #WriteForces true > #WriteMullikenPop 1 > > > > > # 6.11 Options for saving/reading information > UseSaveData true > MD.UseSaveCG true > > #SaveRho true > #SaveDeltaRho true > #SaveElectrostaticPotential true > #SaveTotalPotential true > #SaveIonicCharge true > #SaveTotalCharge true > > %block LocalDensityOfStates > -7.15 12.00 eV > %endblock LocalDensityOfStates > > %block ProjectedDensityOfStates > -7.0 12.0 0.05 521 eV #min, max, broaden 100meV, steps every > 25meV > %endblock ProjectedDensityOfStates > > > #WriteDenchar true > > #%block WaveFuncKPoints > # 0.0 0.0 0.0 from 1040 to 1055 # traps # HOMO 1042 > #%endblock WaveFuncKpoints > > > #OpticalCalculation true > #Optical.EnergyMinimum 0 eV > #ptical.EnergyMaximum 4.5 eV > #Optical.Broaden 0.01 eV > #Optical.NumberOfBands 1100 # LUMO 1085 > > -- ----------------------------------------------------------- Pablo Aguado Puente [email protected] ------------------------------------------------------------
