Thank you Pablo. The problem is this particular pervoskite does not have a centrosymmetric phase. I do realize that polarization is defined modulo x (since momentum operator defined modulo x), then i have to a have reference phase(which i dont). Thanks for the corrections for the PolarizationGrid block. Will try it.
With Best Regards, Andrei Buin. On Tue, Oct 1, 2013 at 9:13 AM, Pablo Aguado <[email protected]> wrote: > > 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] > ------------------------------------------------------------ >
