Is this perovskite ferroelectric? If so you could calculate the
polarization along the switching path (even if it does not go through a
high symmetry phase) to get a meaningful value of P.

Pablo


On Tue, Oct 1, 2013 at 3:20 PM, Andrei Buin <[email protected]> wrote:

> 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]
>> ------------------------------------------------------------
>>
>
>


-- 
-----------------------------------------------------------
Pablo Aguado Puente
[email protected]
------------------------------------------------------------

Responder a