| 
|  Re: open shell systems ( No.1 ) |  |  Date: 2005/12/04 02:57 Name: T.Ozaki
 
Hi, 
 > Question: is this set up for a charged molecule calculation really correct (I seem to
 > add one extra electron to the system, at least in the Atoms.SpeciesAndCoordinates
 > section of the .dat file)?
 
 Your setting is correct. The total number of electons should be neutral
 in the specification of "Atoms.SpeciesAndCoordinates", where each number is
 used for only the generation of the initial charge density. You can control
 the actual total number electons by the keyword, "scf.system.charge".
 Also you can easily check the total number of electrons by summing the Mulliken
 charge up after the calculation to confirm it.
 
 > If the settings are correct, together with set.system.charge = 1, I seem to get wrong
 > dipole momentum of the ammonium ion: it is of the same order as the dipole momentum of
 > the ammonia molecule (nh3). As far as I can grasp, nh4+ is a symmetrical molecule and
 > can not have any sizable dipole momentum (other than from numerical error).
 > The situation with so4(2-) ion seems to be even worse: the dipole momentum is very
 > large, almost 1e3, which is completely unbelievable. Most probably I my setting are
 > too wrong here.
 
 I have also tried the two systems and obtained completely the zero dipole moment
 for both the systems as we expect. Possible source of errors you met can be listed
 as folllows:
 
 (1) Could you obtain the convergent results ?
 (2) scf.lapack.dste=dstevx is better for a symmetric system.
 The lapack routine "dstedc" and "dstegr" tends to fail to correctly calculate
 symmetric systems.
 (3) Remind that OpenMX uses the regular real space grid for the numerical integration.
 So, the total energy and the dipole moment depend on the rotation of system
 if you use a lower cutoff energy (scf.energycutoff), while, of course, the dependency
 on the rotation of system is reduced, as the cutoff energy increases.
 To avoid this dependency, I constructed the structures by modifying "Methane.dat"
 file in the "work" directory, since the relative structural orientation to the regular
 real space grid is symmetrical.
 
 Also, you will see that OpenMX gives a reasonable dipole moment for many molecules,
 and how the dipole moment depends on the choice of basis functions at
 http://netserver.aip.org/cgi-bin/epaps?ID=E-JCPSA6-121-301439
 
 > Question 2: if the settings are right, does it mean that the "default" dft method in
 > nh3.dat file example can handle only closed electronic shells and can not be used
 > for "open shell" calculations, such as molecular ions?
 
 You can try both "closed electronic shells" (spin-unpolarized) and "open shell"
 (spin-polarized) calculations.
 
 Best regards,
 
 T.Ozaki
 
  |  |  Re: open shell systems ( No.2 ) |  |  Date: 2005/12/05 16:35 Name: peter
 
Dear Prof. Ozaki,
 below I am attaching the .dat file for my nh4 calculation. I tried to modify the settings according to your reccomendations, but failed to obtain zero dipole momentum. In fact my result for the dipole momentum (8.12) is very stable with respect to the settings modifications. Please have a look at the .dat file:
 -----------------------------------------------------------
 #
 # SCF calculation of a methane molecule by the LDA
 # and the cluster method
 #
 
 #
 # File Name
 #
 System.CurrrentDirectory ./
 System.Name NH4
 level.of.stdout 1 # default=1 (1-3)
 level.of.fileout 1 # default=1 (0-2)
 
 #
 # Definition of Atomic Species
 #
 
 Species.Number 2
 <Definition.of.Atomic.Species
 H H4.0-s1 H_TM
 N N4.5-s1p1 N_TM_PCC
 Definition.of.Atomic.Species>
 
 #
 # Atoms
 #
 Atoms.Number 5
 Atoms.SpeciesAndCoordinates.Unit Ang # Ang|AU
 <Atoms.SpeciesAndCoordinates
 1 N -0.976418 -0.975766 -0.97872 2.5 2.5 0.0 0.0 1
 2 H -0.965599 0.057111 -0.990461 0.5 0.5 0.0 0.0 1
 3 H -0.00646 -1.32998 -0.950363 0.5 0.5 0.0 0.0 1
 4 H -1.48788 -1.30529 -0.143907 0.5 0.5 0.0 0.0 1
 5 H -1.44574 -1.3249 -1.83015 0.5 0.5 0.0 0.0 1
 Atoms.SpeciesAndCoordinates>
 Atoms.UnitVectors.Unit Ang # Ang|AU
 <Atoms.UnitVectors
 10.0 0.0 0.0
 0.0 10.0 0.0
 0.0 0.0 10.0
 Atoms.UnitVectors>
 
 #
 # SCF or Electronic System
 #
 
 scf.XcType LDA # LDA|LSDA-CA|LSDA-PW|GGA-PBE
 scf.SpinPolarization off # On|Off|NC
 scf.ElectronicTemperature 300.0 # default=300 (K)
 scf.energycutoff 200.0 # default=150 (Ry)
 scf.maxIter 100 # default=40
 scf.EigenvalueSolver cluster # DC|GDC|Cluster|Band
 scf.lapack.dste dstevx # dstegr|dstedc|dstevx, default=dstegr
 scf.Kgrid 1 1 1 # means n1 x n2 x n3
 scf.Mixing.Type rmm-diis # Simple|Rmm-Diis|Gr-Pulay|Kerker|Rmm-Diisk
 scf.Init.Mixing.Weight 0.30 # default=0.30
 scf.Min.Mixing.Weight 0.001 # default=0.001
 scf.Max.Mixing.Weight 0.400 # default=0.40
 scf.Mixing.History 7 # default=5
 scf.Mixing.StartPulay 5 # default=6
 scf.criterion 1.0e-8 # default=1.0e-6 (Hartree)
 
 scf.system.charge 1
 
 #
 # 1D FFT
 #
 
 1DFFT.NumGridK 900 # default=900
 1DFFT.NumGridR 900 # default=900
 1DFFT.EnergyCutoff 2500.0 # default=3DFFT.EnergyCutoff*3.0 (Ry)
 
 #
 # Orbital Optimization
 #
 
 orbitalOpt.Method off # Off|Unrestricted|Restricted
 orbitalOpt.InitCoes Symmetrical # Symmetrical|Free
 orbitalOpt.initPrefactor 0.1 # default=0.1
 orbitalOpt.scf.maxIter 15 # default=12
 orbitalOpt.MD.maxIter 7 # default=5
 orbitalOpt.per.MDIter 20 # default=1000000
 orbitalOpt.criterion 1.0e-4 # default=1.0e-4 (Hartree/borh)^2
 
 #
 # output of contracted orbitals
 #
 
 CntOrb.fileout off # on|off, default=off
 Num.CntOrb.Atoms 1 # default=1
 <Atoms.Cont.Orbitals
 1
 Atoms.Cont.Orbitals>
 
 #
 # SCF Order-N
 #
 
 orderN.HoppingRanges 6.0 # default=5.0 (Ang)
 orderN.NumHoppings 2 # default=2
 
 #
 # MD or Geometry Optimization
 #
 
 MD.Type Opt # Nomd|Opt|NVE|NVT_VS|NVT_NH
 # Constraint_Opt|DIIS
 MD.maxIter 1 # default=1
 MD.TimeStep 1.0 # default=0.5 (fs)
 MD.Opt.criterion 1.0e-4 # default=1.0e-4 (Hartree/bohr)
 
 #
 # restart using a restart file, *.rst
 #
 
 scf.restart off # on|off,default=off
 
 #
 # MO output
 #
 
 MO.fileout off # on|off
 num.HOMOs 1 # default=1
 num.LUMOs 1 # default=1
 MO.Nkpoint 1 # default=1
 <MO.kpoint
 0.0 0.0 0.0
 MO.kpoint>
 
 #
 # DOS and PDOS
 #
 
 Dos.fileout off # on|off, default=off
 Dos.Erange -10.0 10.0 # default = -20 20
 Dos.Kgrid 1 1 1 # default = Kgrid1 Kgrid2 Kgrid3
 
 #
 # output Hamiltonian and overlap
 #
 
 HS.fileout off # on|off, default=off
 
 
  |  |  Re: open shell systems ( No.3 ) |  |  Date: 2005/12/05 16:40 Name: peter
 
Please have a look at the end of the openmx output for the .dat file above:----------------------
 ******************* MD= 1 SCF=11 *******************
 <Poisson> Poisson's equation using FFT...
 <Set_Hamiltonian> Hamiltonian matrix for VNA+dVH+Vxc...
 <Cluster> Eigenvalue problem...
 Atom 1 MulP 2.9671 2.9671 sum 5.9343
 Atom 2 MulP 0.2582 0.2582 sum 0.5163
 Atom 3 MulP 0.2582 0.2582 sum 0.5164
 Atom 4 MulP 0.2582 0.2582 sum 0.5165
 Atom 5 MulP 0.2582 0.2582 sum 0.5165
 <DFT> Total spin S = 0.000000000000
 <DFT> Mixing_weight= 0.348032419010
 <DFT> Uele = -6.086293721744 dUele = 0.000000002256
 <DFT> NormRD = 0.000000000000 Criterion = 0.000000010000
 <MD= 1> Force calculation
 Force calculation #1
 Force calculation #2
 Force calculation #3
 Force calculation #5
 <MD= 1> Double-counting Energy
 Force calculation #6
 Force calculation #7
 
 *******************************************************
 Dipole moment (Debye)
 *******************************************************
 
 Absolute D 8.12725258
 
 Dx Dy Dz
 Total -4.68953265 -4.68622176 -4.70104709
 Core -42.20947180 -42.18122895 -42.30895582
 Electron 37.51993916 37.49500719 37.60790873
 
 *******************************************************
 Kohn-Sham's energies (Hartree) at MD = 1
 *******************************************************
 
 Uele = -6.086293721744
 Uxc0 = 0.400301653350
 Uxc1 = 0.400301653350
 UH0 = -20.588765127127
 UH1 = 1.578292054581
 Ucore = 12.127643937873
 Utot = -12.168519549718
 
 *******************************************************
 Computational times (s) at MD = 1
 *******************************************************
 
 Set_OLP_Kin = 0.17895
 Set_Nonlocal = 0.40981
 Set_Hamiltonian = 0.71516
 Poisson = 2.76630
 diagonalization = 0.00572
 Mixing_DM = 0.00060
 Force = 0.87290
 Correction_Energy = 0.65550
 
 *******************************************************
 MD or geometry opt. at MD = 1
 *******************************************************
 
 <Steepest_Descent> SD_scaling= 0.930797663119
 <Steepest_Descent> |Maximum force| (Hartree/Bohr) = 0.088589990501
 <Steepest_Descent> Criterion (Hartree/Bohr) = 0.000100000000
 
 atom= 1, XYZ(ang) Fxyz(a.u.)= -0.9764 -0.9758 -0.9787 0.0000 0.0000 -0.0000
 atom= 2, XYZ(ang) Fxyz(a.u.)= -0.9656 0.0607 -0.9905 -0.0009 -0.0886 0.0010
 atom= 3, XYZ(ang) Fxyz(a.u.)= -0.0031 -1.3312 -0.9503 -0.0830 0.0303 -0.0024
 atom= 4, XYZ(ang) Fxyz(a.u.)= -1.4897 -1.3064 -0.1410 0.0438 0.0282 -0.0715
 atom= 5, XYZ(ang) Fxyz(a.u.)= -1.4474 -1.3261 -1.8331 0.0402 0.0299 0.0729
 
 outputting data on grids to files...
 ---------------------------------------------
 
 
 As you can see, the molecule has a very large dipole momentum.
 
 Another experiment I did is the following. I took your Methane.dat file from the work directory of your distribution. After scf.system.charge=1 setting, openmx finds the dipole moment 5.4e-5 (vs. 1.5e-7 for the neutral molecule). The results seems to be fine, does it mean that a small change in the orientation of the molecule (Methane.dat vs. our nh4.dat) can lead to such a drastic error in the scf calculation?
 
 Could you please provide me with your .dat file for nh4+ calculation? Are there other reccomendations we could follow?
 
 Thank you very much for your promt answers,
 Peter
  |  |  Re: open shell systems ( No.4 ) |  |  Date: 2005/12/06 01:18 Name: peter
 
note added:
 when nitrogen atom is placed instead of carbon right in the methane.dat file, the dipole momentum of the nh4 molecule appears to be zero after openmx calculation. If a molecular dynamics is used to do geometry optimization, the bond length is fixed, but a small dipole momentum appears.
 
 I tried to calculate hydroxile ion oh-. The .dat file for the calculation i obtained from your water.dat by removing hydrogen atom. The dipole momentum of the ion appears 0.68-0.77 depending on the basis set and a few settings. The "experimental value" can be found on the net and are all larger than one, e.g. 1.58D (http://www.colby.edu/chemistry/webmo/OH-.html).
 
 An attempt to calculate ch3o- ion also leads to a giant dipole momentum, while the dipole momentum for the neutral molecule (ch3oh) is fine.
 
 What is the source for this instability? Does it mean that all ions calculations are unstable with respect to the relative orientation of a molecular ion vs. the FFT "cavity"?
 
 regards,
 Peter
  |  |  Re: open shell systems ( No.5 ) |  |  Date: 2005/12/06 02:35 Name: T.Ozaki
 
Hi,
 Thanks for your detailed reports.
 
 I noticed from your calculations a thing I did not notice
 for calculation of the dipole moment of charged systems.
 
 For a charged system, the dipole moment is not well defined
 if the background contribution is not taken into account.
 
 For instance, imagine a charged linear molecule "A1-A2-A3" along x-axis,
 whose x-coordinate is X-d, X, and X+d for A1, A2, and A3, respectively.
 And the charge is Z, 0, and Z for A1, A2, and A3, respectively.
 Then, the dipole moment is calculated as
 
 Z(X+d) + 0(X) + Z(X-d) = 2ZX
 
 We cleary see that the dipole moment of the charged molecule depends on
 the coordinate X, and in fact this is the source of the problem you met,
 while the background contribution was zero (?un)fortunately in my structures.
 
 The problem is resolved by taking into account the contribution
 of the back ground charge. Then, the total system becomes neutral, and
 the coordinate dependency of the dipole moment disappears.
 
 The following is the result including the contribution of the background
 charge calculated by your input file.
 
 *******************************************************
 Dipole moment (Debye)
 *******************************************************
 
 Absolute D 0.00071255
 
 Dx Dy Dz
 Total 0.00041165 0.00057932 -0.00005157
 Core -42.20947180 -42.18122895 -42.30895582
 Electron 37.51993916 37.49500719 37.60790873
 Back ground 4.68994430 4.68680108 4.70099552
 
 We see that the back ground contribution compensates other contributions.
 Please download the following two modified files, in which the contribution
 of the background charge is taken into account, and replace the corresponding
 files in OpenMX2.73 (note that the two files are valid for only OpenMX2.73).
 http://staff.aist.go.jp/t-ozaki/bugfixed/05Dec06/Correction_Energy.c
 http://staff.aist.go.jp/t-ozaki/bugfixed/05Dec06/openmx_common.h
 
 Regards,
 
 T.Ozaki
 
  |  |  Re: open shell systems ( No.6 ) |  |  Date: 2005/12/08 00:56 Name: peter
 
Dear Prof. Ozaki,
 thank you for your explanations. I played with the settings in the . dat file and managed to find very good dipole momentum for small organic molecules. The next test I did is the molecular polarization. I applied the electric field in the orthogonal directions and derived the polarizability tensor for a few simple molecules (up to 6 atoms). The common observation is that the polarizability seems to be 30-50% lower than I expected from experimental data. This holds even for molecules with good numerical values of dipole momentum. To be specific: I diagonalized the polarizability tensor to compare with the triples of polarizability values given in the literature. If the tensor components where not found, I used the mean value.
 
 Polarizability must be sensitive to inclusion of higher AO in the basis set. Have you performed polarizability calculations with openmx? Is there a place I could read about it and figure out the proper settings?
 
 thank you for your answers,
 Peter
  |  |  Re: openmx runtime questions ( No.7 ) |  |  Date: 2005/12/08 20:01 Name: T.Ozaki
 
Hi,
 The calculation of polarizability has never been reported using OpenMX.
 I guess that a relatively high quality basis set will be required
 for the convergent result.
 
 Regards,
 
 T.Ozaki
 |  |