Siesta

Merge lp:~albertog/siesta/4.1-xc into lp:siesta/4.1

  1. 4.1-xc
  2. Merge into rel-4.1
Proposed by Alberto Garcia
Status: Merged
Merged at revision: 1066
Proposed branch: lp:~albertog/siesta/4.1-xc
Merge into: lp:siesta/4.1
Diff against target: 4146 lines (+351/-3074) (has conflicts)
16 files modified
Docs/developer/ford-pages/grid.md (+1/-1)
Docs/siesta.tex (+14/-6)
Src/Makefile (+26/-30)
Src/SiestaXC/cellxc.F90 (+18/-1)
Src/atom.F (+16/-14)
Src/bsc_cellxc.F (+25/-8)
Src/bsc_xcmod.F (+0/-115)
Src/cellxc_mod.F (+0/-47)
Src/dhscf.F (+143/-145)
Src/forhar.F (+60/-101)
Src/ldau.F (+0/-8)
Src/meshsubs.F (+6/-4)
Src/siesta_init.F (+0/-7)
Src/xc.f (+0/-2547)
Util/Gen-basis/Makefile (+36/-40)
version.info (+6/-0)
Text conflict in version.info
To merge this branch: bzr merge lp:~albertog/siesta/4.1-xc
Reviewer Review Type Date Requested Status
Nick Papior Needs Information
Review via email: mp+361062@code.launchpad.net

Commit message

Simplify the choice between SiestaXC and BSC's versions of cellxc

The BSC_CELLXC preprocessor blocks have been removed and replaced
by conditional blocks which can be selected at runtime.

The user can control which version of cellxc is used by means of the
fdf logical variable:

     XC.Use.BSC.CellXC (T/F)

which is 'false' by default, thus using SiestaXC's cellxc.

The BSC version might be slightly better for GGA operations, and the
SiestaXC version must of course be used when dealing with van der Waals
density functionals.

NOTES:

* BSC's version of cellxc now uses ldaxc, ggaxc, and setxc/getxc from
  SiestaXC. This enhances its functionality (more functionals) and
  saves on duplicated code. The old 'xc.f' file has been removed.

* In 'meshsubs', 'distriphionmesh' now accepts an extra argument to
  flag the need for stencil initilization.

* In 'forhar', a number of 'reord' operations towards 'sequential'
  fine-point ordering (which were compiled-in only for the new
  interface) have been removed, as all the arrays involved are already
  in 'sequential' form.

  Tests show that these changes do not seem to affect the results, however.

  The 'forhar' interface uses the 'ntm' argument in all cases, to simplify
  the code.

* Removed bogus incompatibility check in ldau.F

To post a comment you must log in.
Revision history for this message
Nick Papior (nickpapior) wrote :

Looks good to me!

1) You mention BSC may be "better" for GGA, I guess this means performance wise, no?

2) It is weird that forhar does not show any differences when using the "wrong" ordering? Why, it seems to suggest something is wrong! ? Or?

Lastly, a suggestion on the name of the XC.Use.BSC.CellXC name.

Could we perhaps make this more generalized so one can later use this to select distribution options?

I.e.

   XC.Distribution BSC|GridXC

or something else?

As far as I understand the only difference is the distribution of mesh-points? Correct me if I am wrong.

Great job!

review: Needs Information
Revision history for this message
Nick Papior (nickpapior) wrote :

See additional comments.

lp:~albertog/siesta/4.1-xc updated
1052. By Alberto Garcia

Clarifications in forhar + removal of dead code

1053. By Alberto Garcia

Use "linear" ("yes/no") distribution as seed in SiestaXC

* Give SiestaXC a head start with the appropriate distribution by
  using the "yes/no" (called 'LINEAR') distribution from BSC as the
  initial one. Just specify the box (with the appropriate conversions)
  and use the new cellXC option described below to deactivate further
  internal changes.

* By settting the optional argument 'keep_input_distribution' to
  .true., the internal workload balancer in cellXC is deactivated,
  and the input distribution kept for the XC computations. This is
  advantageous when the client program has better workload-balancers
  available.

1054. By Alberto Garcia

Update grid.md in ford-pages

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'Docs/developer/ford-pages/grid.md'
2--- Docs/developer/ford-pages/grid.md 2018-10-12 23:51:05 +0000
3+++ Docs/developer/ford-pages/grid.md 2019-01-30 13:43:54 +0000
4@@ -57,7 +57,7 @@
5 every point. A ``uniform`` distribution is appropriate. Each process
6 handles a parallelepipedic portion of the mesh box of the same size.
7
8-* The XC routine does not work on points where there is no charge
9+* The XC routine has no work to do on points where there is no charge
10 density. A special algorithm determines the allocation of points to
11 processors so that all have approximately the same workload. This is
12 the ``linear`` distribution (this is a misnomer: it should be
13
14=== modified file 'Docs/siesta.tex'
15--- Docs/siesta.tex 2019-01-30 11:09:09 +0000
16+++ Docs/siesta.tex 2019-01-30 13:43:54 +0000
17@@ -3703,6 +3703,15 @@
18
19 \end{fdfentry}
20
21+\begin{fdflogicalF}{XC!Use.BSC.CellXC}
22+ \index{exchange-correlation!cellXC}
23+
24+ If \fdftrue, the version of \texttt{cellXC} from the BSC's mesh
25+ suite is used instead of the default SiestaXC version. BSC's version
26+ might be slightly better for GGA operations. SiestaXC's version is
27+ mandatory when dealing with van der Waals functionals.
28+
29+\end{fdflogicalF}
30
31
32 \subsection{Spin polarization}
33@@ -12587,12 +12596,11 @@
34 generation of the basis orbitals and the screened pseudopotentials.
35
36 \item%
37- The exchange-correlation routines contained in file xc.f were
38- written by J.M.Soler in 1996 and 1997, in collaboration with
39- \textbf{C.\ Balb\'as} and \textbf{J. L.\ Martins}. Routine pzxc (in
40- the same file), which implements the Perdew-Zunger LDA
41- parametrization of xc, is based on routine velect, written by
42- \textbf{S.\ Froyen}.
43+ The exchange-correlation routines contained in SiestaXC were written
44+ by J.M.Soler in 1996 and 1997, in collaboration with
45+ \textbf{C.\ Balb\'as} and \textbf{J. L.\ Martins}. Routine pzxc,
46+ which implements the Perdew-Zunger LDA parametrization of xc, is
47+ based on routine velect, written by \textbf{S.\ Froyen}.
48
49 \item%
50 The serial version of the multivariate fast fourier transform used
51
52=== modified file 'Src/Makefile'
53--- Src/Makefile 2019-01-23 08:12:54 +0000
54+++ Src/Makefile 2019-01-30 13:43:54 +0000
55@@ -81,7 +81,7 @@
56
57 OBJS = automatic_cell.o atom_options.o \
58 arw.o atomlwf.o bands.o basis_enthalpy.o bessph.o bonds.o \
59- born_charge.o cellxc_mod.o cgwf.o chkdim.o chkgmx.o \
60+ born_charge.o cgwf.o chkdim.o chkgmx.o \
61 chempot.o coceri.o coxmol.o cross.o compute_norm.o\
62 denmat.o denmatlomem.o detover.o dfscf.o diagon.o digcel.o \
63 fft.o dhscf.o constr.o diagk_file.o \
64@@ -150,7 +150,7 @@
65 moreParallelSubs.o \
66 read_xc_info.o \
67 siesta_master.o \
68- bsc_xcmod.o bsc_cellxc.o xc.o \
69+ bsc_cellxc.o \
70 vacuum_level.o \
71 write_orb_indx.o \
72 die.o m_pexsi.o m_pexsi_driver.o m_pexsi_dos.o m_pexsi_local_dos.o\
73@@ -543,8 +543,8 @@
74 atm_transfer.o: periodic_table.o radial.o sys.o
75 atm_types.o: precision.o radial.o
76 atmfuncs.o: atm_types.o precision.o radial.o spher_harm.o sys.o
77-atom.o: atmparams.o atom_options.o basis_specs.o basis_types.o bsc_xcmod.o
78-atom.o: interpolation.o m_filter.o old_atmfuncs.o periodic_table.o precision.o
79+atom.o: atmparams.o atom_options.o basis_specs.o basis_types.o interpolation.o
80+atom.o: m_filter.o old_atmfuncs.o periodic_table.o precision.o
81 atom.o: pseudopotential.o sys.o
82 atom_graph.o: alloc.o atm_types.o atmfuncs.o class_OrbitalDistribution.o
83 atom_graph.o: class_SpData1D.o class_SpData2D.o class_SpData2D.o
84@@ -573,16 +573,14 @@
85 broadcast_projections.o: m_trialorbitalclass.o parallel.o siesta2wannier90.o
86 broyden_optim.o: m_broyddj_nocomm.o m_memory.o parallel.o precision.o
87 broyden_optim.o: siesta_options.o sys.o units.o
88-bsc_cellxc.o: alloc.o bsc_xcmod.o cellxc_mod.o mesh.o moremeshsubs.o parallel.o
89-bsc_cellxc.o: parallelsubs.o precision.o sys.o
90-bsc_xcmod.o: parallel.o precision.o sys.o
91+bsc_cellxc.o: alloc.o mesh.o moremeshsubs.o parallel.o parallelsubs.o
92+bsc_cellxc.o: precision.o sys.o
93 cart2frac.o: precision.o sys.o
94 cell_broyden_optim.o: m_broyddj_nocomm.o parallel.o precision.o sys.o units.o
95 cell_broyden_optim.o: zmatrix.o
96 cell_fire_optim.o: alloc.o m_fire.o m_memory.o parallel.o precision.o
97 cell_fire_optim.o: siesta_options.o sys.o zmatrix.o
98 cellsubs.o: precision.o
99-cellxc_mod.o: bsc_xcmod.o sys.o
100 cgvc.o: alloc.o conjgr_old.o m_mpi_utils.o parallel.o precision.o sys.o units.o
101 cgvc_zmatrix.o: alloc.o conjgr.o m_mpi_utils.o parallel.o precision.o sys.o
102 cgvc_zmatrix.o: units.o zmatrix.o
103@@ -672,13 +670,13 @@
104 detover.o: alloc.o linpack.o parallel.o parallelsubs.o precision.o
105 dfscf.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o mesh.o meshdscf.o
106 dfscf.o: meshphi.o parallel.o parallelsubs.o precision.o sys.o
107-dhscf.o: alloc.o atmfuncs.o bsc_xcmod.o cellxc_mod.o delk.o dfscf.o
108-dhscf.o: doping_uniform.o files.o forhar.o iogrid_netcdf.o m_charge_add.o
109-dhscf.o: m_efield.o m_hartree_add.o m_iorho.o m_mesh_node.o m_ncdf_io.o
110-dhscf.o: m_ncdf_siesta.o m_partial_charges.o m_rhog.o m_spin.o
111-dhscf.o: m_ts_global_vars.o m_ts_hartree.o m_ts_options.o m_ts_voltage.o mesh.o
112-dhscf.o: meshdscf.o meshsubs.o moremeshsubs.o parallel.o parsing.o precision.o
113-dhscf.o: rhofft.o rhoofd.o siesta_options.o sys.o units.o vmat.o
114+dhscf.o: alloc.o atmfuncs.o delk.o dfscf.o doping_uniform.o files.o forhar.o
115+dhscf.o: iogrid_netcdf.o m_charge_add.o m_efield.o m_hartree_add.o m_iorho.o
116+dhscf.o: m_mesh_node.o m_ncdf_io.o m_ncdf_siesta.o m_partial_charges.o m_rhog.o
117+dhscf.o: m_spin.o m_ts_global_vars.o m_ts_hartree.o m_ts_options.o
118+dhscf.o: m_ts_voltage.o mesh.o meshdscf.o meshsubs.o moremeshsubs.o parallel.o
119+dhscf.o: parsing.o precision.o rhofft.o rhoofd.o siesta_options.o sys.o units.o
120+dhscf.o: vmat.o
121 diag.o: alloc.o diag_option.o parallel.o precision.o sys.o
122 diag2g.o: fermid.o intrinsic_missing.o parallel.o parallelsubs.o precision.o
123 diag2g.o: sys.o
124@@ -1104,9 +1102,9 @@
125 meshdscf.o: alloc.o atomlist.o m_dscfcomm.o meshphi.o parallel.o parallelsubs.o
126 meshdscf.o: precision.o
127 meshphi.o: alloc.o precision.o
128-meshsubs.o: alloc.o atm_types.o atmfuncs.o cellsubs.o cellxc_mod.o chkgmx.o
129-meshsubs.o: fft1d.o mesh.o meshphi.o moremeshsubs.o parallel.o parallelsubs.o
130-meshsubs.o: precision.o radial.o siesta_cml.o sys.o
131+meshsubs.o: alloc.o atm_types.o atmfuncs.o cellsubs.o chkgmx.o fft1d.o mesh.o
132+meshsubs.o: meshphi.o moremeshsubs.o parallel.o parallelsubs.o precision.o
133+meshsubs.o: radial.o siesta_cml.o sys.o
134 metaforce.o: alloc.o parallel.o precision.o sys.o
135 minvec.o: cellsubs.o precision.o sorting.o sys.o
136 mixer.o: atomlist.o m_mixing.o m_mixing_scf.o m_spin.o parallel.o precision.o
137@@ -1265,17 +1263,16 @@
138 siesta_forces.o: state_analysis.o state_init.o sys.o timer.o units.o
139 siesta_forces.o: write_subs.o
140 siesta_geom.o: precision.o
141-siesta_init.o: alloc.o atomlist.o bands.o bsc_xcmod.o
142-siesta_init.o: class_Fstack_Pair_Geometry_SpData2D.o diag_option.o files.o
143-siesta_init.o: flook_siesta.o ioxv.o kpoint_grid.o kpoint_pdos.o ksvinit.o
144-siesta_init.o: m_check_walltime.o m_cite.o m_energies.o m_eo.o m_fixed.o
145-siesta_init.o: m_forces.o m_iostruct.o m_mpi_utils.o m_new_dm.o m_rmaxh.o
146-siesta_init.o: m_spin.o m_steps.o m_supercell.o m_timer.o m_wallclock.o
147-siesta_init.o: metaforce.o molecularmechanics.o object_debug.o parallel.o
148-siesta_init.o: parallelsubs.o projected_DOS.o siesta_cmlsubs.o siesta_dicts.o
149-siesta_init.o: siesta_geom.o siesta_options.o sparse_matrices.o struct_init.o
150-siesta_init.o: sys.o timer.o timestamp.o ts_init.o units.o writewave.o
151-siesta_init.o: zmatrix.o
152+siesta_init.o: alloc.o atomlist.o bands.o class_Fstack_Pair_Geometry_SpData2D.o
153+siesta_init.o: diag_option.o files.o flook_siesta.o ioxv.o kpoint_grid.o
154+siesta_init.o: kpoint_pdos.o ksvinit.o m_check_walltime.o m_cite.o m_energies.o
155+siesta_init.o: m_eo.o m_fixed.o m_forces.o m_iostruct.o m_mpi_utils.o
156+siesta_init.o: m_new_dm.o m_rmaxh.o m_spin.o m_steps.o m_supercell.o m_timer.o
157+siesta_init.o: m_wallclock.o metaforce.o molecularmechanics.o object_debug.o
158+siesta_init.o: parallel.o parallelsubs.o projected_DOS.o siesta_cmlsubs.o
159+siesta_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
160+siesta_init.o: struct_init.o sys.o timer.o timestamp.o ts_init.o units.o
161+siesta_init.o: writewave.o zmatrix.o
162 siesta_master.o: iopipes.o iosockets.o precision.o sys.o
163 siesta_move.o: atomlist.o broyden_optim.o cell_broyden_optim.o
164 siesta_move.o: cell_fire_optim.o dynamics.o fire_optim.o flook_siesta.o ioxv.o
165@@ -1346,7 +1343,6 @@
166 writewave.o: alloc.o atmfuncs.o atomlist.o densematrix.o diag.o diag_option.o
167 writewave.o: get_kpoints_scale.o kpoint_grid.o m_spin.o parallel.o
168 writewave.o: parallelsubs.o precision.o siesta_geom.o sys.o units.o
169-xc.o: alloc.o bsc_xcmod.o precision.o sys.o
170 xml.o: precision.o
171 zm_broyden_optim.o: m_broyddj_nocomm.o parallel.o precision.o sys.o units.o
172 zm_broyden_optim.o: zmatrix.o
173
174=== modified file 'Src/SiestaXC/cellxc.F90'
175--- Src/SiestaXC/cellxc.F90 2017-09-04 11:10:16 +0000
176+++ Src/SiestaXC/cellxc.F90 2019-01-30 13:43:54 +0000
177@@ -232,7 +232,8 @@
178 CONTAINS ! nothing else but public routine cellXC
179
180 SUBROUTINE cellXC( irel, cell, nMesh, lb1, ub1, lb2, ub2, lb3, ub3, &
181- nSpin, dens, Ex, Ec, Dx, Dc, stress, Vxc, dVxcdD )
182+ nSpin, dens, Ex, Ec, Dx, Dc, stress, Vxc, dVxcdD, &
183+ keep_input_distribution )
184
185 ! Module routines
186 use mesh3D, only: addMeshData ! Accumulates a mesh array
187@@ -303,6 +304,7 @@
188 ! (spin) xc potential
189 real(gp),intent(out),optional:: & ! dVxc/dDens (LDA only)
190 dVxcdD(0:ub1-lb1,0:ub2-lb2,0:ub3-lb3,1:nSpin**2)
191+ logical, intent(in),optional:: keep_input_distribution
192
193 ! Fix the order of the numerical derivatives
194 ! nn is the number of points used in each coordinate and direction,
195@@ -404,6 +406,8 @@
196 type(allocDefaults):: &
197 prevAllocDefaults
198
199+ logical :: keep_input_distr
200+
201 #ifdef DEBUG_XC
202 ! Variables for debugging
203 real(dp):: GDtot(3), q, dqdD, dqdGD(3)
204@@ -415,6 +419,11 @@
205 ! Start time counter
206 call timer_start( myName )
207
208+ keep_input_distr = .false.
209+ if (present(keep_input_distribution)) then
210+ keep_input_distr = keep_input_distribution
211+ endif
212+
213 #ifdef DEBUG_XC
214 ! Initialize udebug variable
215 call setDebugOutputUnit()
216@@ -469,6 +478,12 @@
217 ! Get ID of the I/O distribution of mesh points
218 call setMeshDistr( ioDistr, nMesh, ioBox )
219
220+ if (keep_input_distr) then
221+
222+ myDistr = ioDistr
223+
224+ else
225+
226 ! If nMesh has changed, use input distribution also initially as myDistr
227 if (any(nMesh/=oldMesh)) call setMeshDistr( myDistr, nMesh, ioBox )
228 oldMesh = nMesh
229@@ -519,6 +534,7 @@
230 #endif /* DEBUG_XC */
231 end if ! (nodes>1 .and. timeDisp/timeAvge>maxUnbalance)
232
233+ endif
234 #ifdef DEBUG_XC
235 ! Keep input distribution for the time being
236 ! myDistr = ioDistr
237@@ -535,6 +551,7 @@
238 ! - The parallel mesh distribution is changed internally (even with LDA)
239 ! - We need finite differences (GGA) and have a distributed mesh
240 if (.not.sameMeshDistr(myDistr,ioDistr) .or. (GGA .and. myDistr/=0)) then
241+
242 m11=myBox(1,1); m12=myBox(1,2); m13=myBox(1,3)
243 m21=myBox(2,1); m22=myBox(2,2); m23=myBox(2,3)
244 ns = nSpin ! Just a shorter name
245
246=== modified file 'Src/atom.F'
247--- Src/atom.F 2018-07-13 18:31:43 +0000
248+++ Src/atom.F 2019-01-30 13:43:54 +0000
249@@ -1,3 +1,5 @@
250+
251+
252 !
253 ! Copyright (C) 1996-2016 The SIESTA group
254 ! This file is distributed under the terms of the
255@@ -39,10 +41,10 @@
256 use basis_types, only: basis_def_t
257 use fdf
258 use pseudopotential, only: pseudopotential_t
259-#ifndef BSC_CELLXC
260+
261 use siestaXC, only: atomXC ! XC for a spherical charge
262 use siestaXC, only: getXC ! Returns the XC functional to be used
263-#endif /* ! BSC_CELLXC */
264+
265
266 implicit none
267
268@@ -1975,20 +1977,20 @@
269 C Checking the functional used for exchange-correlation energy.
270 C Written by D. Sanchez-Portal, Aug. 1998
271
272-#ifdef BSC_CELLXC
273- use bsc_xcmod, only : nXCfunc, XCfunc, XCauth
274-
275-#endif /* BSC_CELLXC */
276+
277+
278+
279+
280 C Passed variable
281 character(len=2), intent(in) :: icorr
282
283 C Local variables
284-#ifndef BSC_CELLXC
285+
286 integer :: nf, nXCfunc
287 character(len=20) :: XCauth(10), XCfunc(10)
288-#else /* BSC_CELLXC */
289- integer :: nf
290-#endif /* BSC_CELLXC */
291+
292+
293+
294 character(len=40) :: ps_string
295
296 ps_string = "Unknown atomic XC code"
297@@ -2001,7 +2003,7 @@
298 if (icorr .eq. "wc") ps_string ="GGA WC"
299 if (icorr .eq. "bl") ps_string ="GGA BLYP"
300 if (icorr .eq. "ps") ps_string ="GGA PBEsol"
301-#ifndef BSC_CELLXC
302+
303 if (icorr .eq. "jo") ps_string ="GGA PBEJsJrLO"
304 if (icorr .eq. "jh") ps_string ="GGA PBEJsJrHEG"
305 if (icorr .eq. "go") ps_string ="GGA PBEGcGxLO"
306@@ -2017,7 +2019,7 @@
307
308 C Get functional(s) being used
309 call getXC( nXCfunc, XCfunc, XCauth )
310-#endif /* ! BSC_CELLXC */
311+
312
313 C Loop over functionals
314 do nf = 1,nXCfunc
315@@ -2095,7 +2097,7 @@
316 . 'xc_check: WARNING: Pseudopotential generated with',
317 . trim(ps_string), " functional"
318
319-#ifndef BSC_CELLXC
320+
321 elseif((XCauth(nf).eq.'PW91').and.(XCfunc(nf).eq.'GGA')) then
322
323 write(6,'(a)')
324@@ -2215,7 +2217,7 @@
325 . 'xc_check: WARNING: Pseudopotential generated with',
326 . trim(ps_string), " functional"
327
328-#endif /* ! BSC_CELLXC */
329+
330 else
331
332 write(6,'(a)')
333
334=== modified file 'Src/bsc_cellxc.F'
335--- Src/bsc_cellxc.F 2016-01-25 16:00:16 +0000
336+++ Src/bsc_cellxc.F 2019-01-30 13:43:54 +0000
337@@ -115,12 +115,11 @@
338 C *******************************************************************
339
340 use precision, only : dp, grid_p
341- use bsc_xcmod, only : nXCfunc, XCfunc, XCauth
342- use bsc_xcmod, only : XCweightX, XCweightC
343 use sys, only : die
344 use alloc, only : re_alloc, de_alloc
345- use cellxc_mod, only : GGA
346 use mesh, only : nsm, meshLim
347+ use siestaxc, only : ggaxc, ldaxc
348+ use siestaxc, only : getXC
349 #ifdef MPI
350 C Modules
351 use mesh, only : nmeshg
352@@ -179,7 +178,7 @@
353 #endif
354
355 C Local variables and arrays
356- logical GGAfunctl
357+ logical GGA, GGAfunctl
358 integer I1, I2, I3, IC, IN, IP, IS, IX,
359 & J1, J2, J3, JN, JP(3,-NN:NN), JS, JX,
360 & KS, NPG, nf
361@@ -190,12 +189,17 @@
362 & DVCDN(mspin*mspin), DVXDN(mspin*mspin),
363 & EPSC, EPSX, F1, F2, GD(3,mspin),
364 & VOLCEL, volume, stress(3,3)
365- external GGAXC, LDAXC, RECLAT, VOLCEL
366+ external RECLAT, VOLCEL
367+
368 integer BS(3), iDistr
369 real(grid_p), pointer :: bdensX(:,:,:), bdensY(:,:,:),
370 & bdensZ(:,:,:), bvxcX(:,:,:),
371 & bvxcY(:,:,:), bvxcZ(:,:,:)
372
373+ integer :: nXCfunc
374+ character(len=20) :: XCauth(10), XCfunc(10)
375+ real(dp) :: XCweightX(10), XCweightC(10)
376+
377 #ifdef DEBUG
378 call write_debug( ' PRE cellxc' )
379 #endif
380@@ -207,9 +211,21 @@
381 #endif
382 call timer( 'cellXC', 1 )
383
384-C Check ider
385- if (ider.ne.0 .and. GGA)
386- $ call die('cellXC: ider=1 available only for LDA')
387+ call getXC( nXCfunc, XCfunc, XCauth, XCweightX, XCweightC)
388+
389+ GGA = .false.
390+ do nf = 1,nXCfunc
391+ if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga' ) then
392+ GGA = .true.
393+ endif
394+ if ( XCfunc(nf).eq.'VDW' .or. XCfunc(nf).eq.'vdw') then
395+ call die("BSC's cellxc cannot handle VDW functionals")
396+ endif
397+ enddo
398+
399+ if (ider.ne.0 .and. GGA) then
400+ call die('BSC cellXC: ider=1 available only for LDA')
401+ endif
402
403 C Check value of mspin
404 if (mspin.lt.nspin) then
405@@ -584,6 +600,7 @@
406 #endif
407 endif
408
409+
410 C Loop over all functionals
411 do nf = 1,nXCfunc
412
413
414=== removed file 'Src/bsc_xcmod.F'
415--- Src/bsc_xcmod.F 2016-01-25 16:00:16 +0000
416+++ Src/bsc_xcmod.F 1970-01-01 00:00:00 +0000
417@@ -1,115 +0,0 @@
418-!
419-! Copyright (C) 1996-2016 The SIESTA group
420-! This file is distributed under the terms of the
421-! GNU General Public License: see COPYING in the top directory
422-! or http://www.gnu.org/copyleft/gpl.txt.
423-! See Docs/Contributors.txt for a list of contributors.
424-!
425- module bsc_xcmod
426-
427- use precision, only : dp
428-
429- implicit none
430-
431- integer, parameter :: MaxFunc = 10
432- integer, save :: nXCfunc
433- character(len=10), save :: XCauth(MaxFunc)
434- character(len=10), save :: XCfunc(MaxFunc)
435- real(dp), save :: XCweightX(MaxFunc)
436- real(dp), save :: XCweightC(MaxFunc)
437-
438- public :: nXCfunc, XCauth, XCfunc
439- public :: XCweightX, XCweightC
440- public :: setXC
441- private
442-
443- contains
444-
445- subroutine setXC()
446-C
447-C This subroutine reads the exchange-correlation functional information
448-C
449- use fdf
450- use parallel, only : Node
451- use sys, only : die
452-
453- implicit none
454-
455-C Local variables
456- integer :: n
457- integer :: ni
458- integer :: nn
459- integer :: nr
460-
461- type(block_fdf) :: bfdf
462- type(parsed_line), pointer :: pline
463-
464-C Read XC functionals
465- if (fdf_block('xc.hybrid',bfdf)) then
466- if (.not. fdf_bline(bfdf,pline)) then
467- call die('setXC: ERROR no data in XC.hybrid block')
468- endif
469- ni = fdf_bnintegers(pline)
470-
471- if (ni .eq. 0) then
472- call die('setXC: Number of functionals missing in ' //
473- . 'XC.hybrid')
474- endif
475- nXCfunc = abs(fdf_bintegers(pline,1))
476- if (nXCfunc .gt. MaxFunc) then
477- call die('setXC: Too many functionals in XC.hybrid')
478- endif
479- do n= 1, nXCfunc
480- if (.not. fdf_bline(bfdf,pline)) then
481- call die('setXC: Number of XC functionals does not match')
482- endif
483- nn = fdf_bnnames(pline)
484- nr = fdf_bnreals(pline)
485-
486- if (nn .gt. 0) then
487- XCfunc(n) = fdf_bnames(pline,1)
488- else
489- XCfunc(n) = 'LDA'
490- endif
491- if (nn .gt. 1) then
492- XCauth(n) = fdf_bnames(pline,2)
493- else
494- XCauth(n) = 'PZ'
495- endif
496- if (nr .gt. 1) then
497- XCweightX(n) = fdf_breals(pline,1)
498- XCweightC(n) = fdf_breals(pline,2)
499- elseif (nr .eq. 1) then
500- XCweightX(n) = fdf_breals(pline,1)
501- XCweightC(n) = fdf_breals(pline,1)
502- else
503- XCweightX(n) = 1.0_dp
504- XCweightC(n) = 1.0_dp
505- endif
506- enddo
507- else
508- nXCfunc = 1
509- XCfunc(1) = fdf_string('xc.functional','LDA')
510- XCauth(1) = fdf_string('xc.authors','PZ')
511- XCweightX(1) = 1.0_dp
512- XCweightC(1) = 1.0_dp
513- endif
514-
515-C Output data for hybrid functionals
516- if ((nXCfunc .gt. 1) .and. (Node .eq. 0)) then
517- write(6,'(/,''xc:'')')
518- write(6,'(''xc: Hybrid exchange-correlation functional:'')
519- . ')
520- write(6,'(''xc:'')')
521- write(6,'(''xc: Number Functional Authors '',
522- . '' Weight(Ex) Weight(Ec)'')')
523- do n = 1,nXCfunc
524- write(6,'(''xc: '',i4,7x,a3,5x,a10,3x,f5.3,8x,f5.3)')
525- . n,XCfunc(n),XCauth(n),XCweightX(n),XCweightC(n)
526- enddo
527- write(6,'(''xc:'')')
528- endif
529-
530- end subroutine setXC
531-
532- end module bsc_xcmod
533
534=== removed file 'Src/cellxc_mod.F'
535--- Src/cellxc_mod.F 2016-01-25 16:00:16 +0000
536+++ Src/cellxc_mod.F 1970-01-01 00:00:00 +0000
537@@ -1,47 +0,0 @@
538-! ---
539-! Copyright (C) 1996-2016 The SIESTA group
540-! This file is distributed under the terms of the
541-! GNU General Public License: see COPYING in the top directory
542-! or http://www.gnu.org/copyleft/gpl.txt .
543-! See Docs/Contributors.txt for a list of contributors.
544-! ---
545- MODULE cellxc_mod
546- PRIVATE
547- PUBLIC :: setGGA, GGA
548- logical :: GGA=.FALSE.
549- CONTAINS
550- subroutine setGGA( )
551-#ifndef BSC_CELLXC
552- use siestaXC, only : getXC ! Returns the XC functional used
553-#else /* BSC_CELLXC */
554- use bsc_xcmod, only : nXCfunc, XCfunc
555-#endif /* BSC_CELLXC */
556- use sys, only : die
557- implicit none
558-C Local variables
559-#ifndef BSC_CELLXC
560- integer :: nf, nXCfunc
561- character(len=20):: XCauth(10), XCfunc(10)
562-#else /* BSC_CELLXC */
563- integer :: nf
564-#endif /* BSC_CELLXC */
565-
566-!---------------------------------------------------------------------- BEGIN
567-#ifndef BSC_CELLXC
568- call getXC( nXCfunc, XCfunc, XCauth )
569-#endif /* ! BSC_CELLXC */
570- do nf = 1,nXCfunc
571- if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then
572- GGA = .true.
573- else
574- if ( XCfunc(nf).ne.'LDA' .and. XCfunc(nf).ne.'lda' .and.
575- & XCfunc(nf).ne.'LSD' .and. XCfunc(nf).ne.'lsd' ) then
576- write(6,*) 'cellXC: Unknown functional ', XCfunc(nf)
577- call die()
578- endif
579- endif
580- enddo
581-!----------------------------------------------------------------------- END
582- end subroutine setGGA
583-
584- END MODULE cellxc_mod
585
586=== modified file 'Src/dhscf.F'
587--- Src/dhscf.F 2018-12-23 19:28:49 +0000
588+++ Src/dhscf.F 2019-01-30 13:43:54 +0000
589@@ -45,6 +45,8 @@
590 & dvol, field(3), rmax, scell(3,3)
591 real(dp) :: G2mesh = 0.0_dp
592 logical :: debug_dhscf = .false.
593+ logical :: use_bsc_cellxc
594+
595 character(len=*), parameter :: debug_fmt =
596 & '('' -- dhscf : Node '',i3,tr1,a,4(tr1,e20.7))'
597
598@@ -66,11 +68,7 @@
599 use sys, only : die
600 use mesh, only : xdsp, nsm, nsp, meshLim
601 use parsing
602-#ifndef BSC_CELLXC
603- use siestaXC, only : getXC ! Returns the XC functional used
604-#else /* BSC_CELLXC */
605- use bsc_xcmod, only : nXCfunc, XCauth
606-#endif /* BSC_CELLXC */
607+
608 use alloc, only : re_alloc, de_alloc
609 use siesta_options, only : harrisfun
610 use meshsubs, only : PartialCoreOnMesh
611@@ -93,9 +91,6 @@
612 use m_ncdf_siesta, only : cdf_init_grid
613 use m_ncdf_io, only : cdf_init_mesh
614 #endif
615-#ifdef BSC_CELLXC
616- use cellxc_mod, only : setGGA
617-#endif /* BSC_CELLXC */
618 use m_efield, only : initialize_efield, acting_efield
619 use m_efield, only : get_field_from_dipole
620 use m_efield, only : dipole_correction
621@@ -145,10 +140,8 @@
622 integer, pointer :: numphi(:), numphi_par(:)
623
624 integer :: nm(3) ! For call to initMesh
625-#ifndef BSC_CELLXC
626- integer :: nXCfunc
627- character(len=20) :: XCauth(10), XCfunc(10)
628-#endif /* ! BSC_CELLXC */
629+ logical :: need_gradients_in_xc, is_vdw
630+
631 ! Transport direction (unit-cell aligned)
632 integer :: iE
633 real(dp) :: ortho, field(3), field2(3)
634@@ -182,18 +175,6 @@
635 & call die('dhscf: ERROR: spiral defined but nspin < 4')
636 endif ! First time
637
638-#ifndef BSC_CELLXC
639-C Get functional(s) being used
640- call getXC( nXCfunc, XCfunc, XCauth )
641-#endif /* ! BSC_CELLXC */
642-
643- if (harrisfun) then
644- do n = 1,nXCfunc
645- if (.not.(leqi(XCauth(n),'PZ').or.leqi(XCauth(n),'CA'))) then
646- call die("** Harris forces not implemented for non-LDA XC")
647- endif
648- enddo
649- endif
650
651 C ----------------------------------------------------------------------
652 C Orbital initialisation : part 1
653@@ -278,11 +259,23 @@
654 C Initialise quantities relating to the atom-mesh positioning
655 call InitAtomMesh( UNIFORM, na, xa )
656
657-#ifdef BSC_CELLXC
658-C Check if we need extencils in cellxc
659- call setGGA( )
660-
661-#endif /* BSC_CELLXC */
662+! Check if we need stencils in cellxc, and detect incompatibility
663+! with Harris forces
664+ call xc_setup(need_gradients_in_xc, is_vdw)
665+
666+ if (need_gradients_in_xc .and. harrisfun) then
667+ call die("** Harris forces not implemented for non-LDA XC")
668+ endif
669+
670+ ! Give the opportunity to use BSC's version,
671+ ! unless we have a vdw functional
672+ use_bsc_cellxc = fdf_get("XC.Use.BSC.CellXC", .false.)
673+ if (is_vdw .and. use_bsc_cellxc) then
674+ call message("INFO",
675+ $ "BSC's cellxc cannot be used with vdW functionals")
676+ use_bsc_cellxc = .false.
677+ endif
678+
679 C Compute the number of orbitals on the mesh and recompute the
680 C partions for every processor in order to have a similar load
681 C in each of them.
682@@ -295,7 +288,7 @@
683 !$OMP end parallel do
684
685 call distriPhiOnMesh( nm, nmpl, norb, iaorb, iphorb,
686- & isa, numphi )
687+ & isa, numphi, need_gradients_in_xc )
688
689 C Find if there are partial-core-corrections for any atom
690 npcc = 0
691@@ -623,18 +616,13 @@
692 C ----------------------------------------------------------------------
693 C Modules
694 use precision, only : dp, grid_p
695-#ifndef BSC_CELLXC
696- use parallel, only : ProcessorY
697-#endif /* ! BSC_CELLXC */
698
699-! Number of Mesh divisions of each cell vector (global)
700-! The status of this variable is confusing
701 use parallel, only : Node, Nodes
702 use atmfuncs, only : rcut, rcore
703 use units, only : Debye, eV, Ang
704 use fdf
705 use sys, only : die, bye
706- use mesh, only : nsm, nsp
707+ use mesh, only : nsm, nsp, meshLim
708 use parsing
709 use m_iorho, only : write_rho
710 use m_forhar, only : forhar
711@@ -652,14 +640,12 @@
712 use moreMeshSubs, only : TO_SEQUENTIAL, TO_CLUSTER, KEEP
713 use m_partial_charges, only: compute_partial_charges
714 use m_partial_charges, only: want_partial_charges
715-#ifndef BSC_CELLXC
716
717- use siestaXC, only : cellXC ! Finds xc energy and potential
718+ use siestaXC, only : cellXC ! Finds xc energy and potential
719 use siestaXC, only : myMeshBox ! Returns my processor mesh box
720 use siestaXC, only : jms_setMeshDistr => setMeshDistr
721 ! Sets a distribution of mesh
722 ! points over parallel processors
723-#endif /* BSC_CELLXC */
724 use m_vmat, only : vmat
725 use m_rhoofd, only: rhoofd
726 #ifdef MPI
727@@ -754,10 +740,9 @@
728 C logical IsDiag : Is supercell diagonal?
729 C integer ispin : Spin index
730 C integer j : General-purpose index
731-#ifndef BSC_CELLXC
732-C integer JDGdistr : J.D.Gale's parallel distribution of mesh points
733-C integer myBox(2,3) : My processor's mesh box
734-#endif /* ! BSC_CELLXC */
735+
736+C integer myBox(2,3) : My processor's mesh box
737+
738 C integer nbcell : Number of independent bulk lattice vectors
739 C integer npcc : Partial core corrections? (0=no, 1=yes)
740 C integer nsd : Number of diagonal spin values (1 or 2)
741@@ -780,22 +765,11 @@
742
743 C Local variables
744 integer :: i, ia, ip, ispin, nsd, np_vac
745-#ifndef BSC_CELLXC
746-! Interface to JMS's SiestaXC
747- integer :: myBox(2,3)
748- integer, save :: JDGdistr=-1
749- real(dp) :: stressXC(3,3)
750-#endif /* ! BSC_CELLXC */
751 real(dp) :: b1Xb2(3), const, DEc, DEx, DStres(3,3),
752 & Ec, Ex, rhotot, x0(3), volume, Vmax_vac, Vmean_vac
753-#ifdef BSC_CELLXC
754-! Dummy arrays for cellxc call
755- real(grid_p) :: aux3(3,1)
756- real(grid_p) :: dummy_DVxcdn(1,1,1)
757-#endif /* BSC_CELLXC */
758
759 logical :: use_rhog
760-
761+
762 real(dp), external :: volcel, ddot
763
764 external
765@@ -805,9 +779,16 @@
766 & reord, rhooda, rhoofdsp,
767 & timer, vmatsp,
768 & readsp
769-#ifdef BSC_CELLXC
770+
771+! Dummy arrays for bsc_cellxc call
772+ real(grid_p) :: aux3(3,1)
773+ real(grid_p) :: dummy_DVxcdn(1,1,1)
774+ real(grid_p), pointer :: Vscf_gga(:,:), DRho_gga(:,:)
775 external bsc_cellxc
776-#endif /* BSC_CELLXC */
777+
778+! Interface to JMS's SiestaXC
779+ integer :: myBox(2,3)
780+ real(dp) :: stressXC(3,3)
781
782 C Work arrays
783 real(grid_p), pointer :: Vscf(:,:), Vscf_par(:,:),
784@@ -818,9 +799,7 @@
785 & DRho_quad(:,:) => null()
786 ! Temporary reciprocal spin quantity
787 real(grid_p) :: rnsd
788-#ifdef BSC_CELLXC
789- real(grid_p), pointer :: Vscf_gga(:,:), DRho_gga(:,:)
790-#endif /* BSC_CELLXC */
791+
792 #ifdef MPI
793 integer :: MPIerror
794 real(dp) :: sbuffer(7), rbuffer(7)
795@@ -862,9 +841,7 @@
796
797 nullify( Vscf, Vscf_par, DRho, DRho_par,
798 & Vaux, Vaux_par, Chlocal, Totchar )
799-#ifdef BSC_CELLXC
800 nullify( Vscf_gga, DRho_gga)
801-#endif /* BSC_CELLXC */
802
803 volume = volcel(cell)
804
805@@ -1471,19 +1448,6 @@
806 !$OMP end parallel do
807 endif
808
809-C ----------------------------------------------------------------------
810-#ifndef BSC_CELLXC
811-C Set uniform distribution of mesh points and find my processor mesh box
812-C This is the interface to JM Soler's own cellxc routine, which sets
813-C up the right distribution internally.
814-C ----------------------------------------------------------------------
815-
816- call jms_setMeshDistr( distrID=JDGdistr, nMesh=ntm,
817- . nNodesX=1, nNodesY=ProcessorY, nBlock=nsm )
818- call myMeshBox( ntm, JDGdistr, myBox )
819-
820-C ----------------------------------------------------------------------
821-#endif /* ! BSC_CELLXC */
822 C Exchange-correlation energy
823 C ----------------------------------------------------------------------
824 call re_alloc( Vscf, 1, ntpl, 1, nspin, 'Vscf', 'dhscf' )
825@@ -1534,64 +1498,83 @@
826
827 ! Everything now is in UNIFORM, sequential form
828
829- call timer("CellXC",1)
830-#ifdef BSC_CELLXC
831- if (nodes.gt.1) then
832- call setMeshDistr( LINEAR, nsm, nsp, nml, nmpl, ntml, ntpl )
833- endif
834-
835- call re_alloc( Vscf_gga, 1, ntpl, 1, nspin, 'Vscf_gga', 'dhscf' )
836- call re_alloc( DRho_gga, 1, ntpl, 1, nspin, 'DRho_gga', 'dhscf' )
837-
838- ! Redistribute all spin densities
839- do ispin = 1, nspin
840- fsrc => DRho(:,ispin)
841- fdst => DRho_gga(:,ispin)
842- call distMeshData( UNIFORM, fsrc, LINEAR, fdst, KEEP )
843- enddo
844-
845- call bsc_cellxc( 0, 0, cell, ntml, ntml, ntpl, 0, aux3, nspin,
846- & DRho_gga, Ex, Ec, DEx, DEc, Vscf_gga,
847- & dummy_DVxcdn, stressl )
848-#endif /* BSC_CELLXC */
849-
850-#ifndef BSC_CELLXC
851- call cellXC( 0, cell, ntm, myBox(1,1), myBox(2,1),
852- . myBox(1,2), myBox(2,2),
853- . myBox(1,3), myBox(2,3), nspin,
854- . DRho, Ex, Ec, DEx, DEc, stressXC, Vscf )
855-#else /* BSC_CELLXC */
856-
857- if (nodes.gt.1) then
858- call setMeshDistr( UNIFORM, nsm, nsp, nml, nmpl, ntml, ntpl )
859- endif
860-
861- ! Redistribute to the Vxc array
862- do ispin = 1,nspin
863- fsrc => Vscf_gga(:,ispin)
864- fdst => Vscf(:,ispin)
865- call distMeshData( LINEAR, fsrc, UNIFORM, fdst, KEEP )
866- enddo
867-#endif /* BSC_CELLXC */
868+ call timer("XC",1)
869+
870+ ! Switch to "zero/not-zero rho" distribution (miscalled 'LINEAR')
871+ if (nodes.gt.1) then
872+ call setMeshDistr( LINEAR, nsm, nsp, nml, nmpl, ntml, ntpl )
873+ endif
874+
875+ call re_alloc( Vscf_gga, 1, ntpl, 1, nspin, 'Vscf_gga','dhscf')
876+ call re_alloc( DRho_gga, 1, ntpl, 1, nspin, 'DRho_gga','dhscf')
877+
878+ ! Redistribute all spin densities
879+ do ispin = 1, nspin
880+ fsrc => DRho(:,ispin)
881+ fdst => DRho_gga(:,ispin)
882+ call distMeshData( UNIFORM, fsrc, LINEAR, fdst, KEEP )
883+ enddo
884+
885+ if (use_bsc_cellxc) then
886+
887+ call timer("BSC-CellXC",1)
888+ call bsc_cellxc( 0, 0, cell, ntml, ntml, ntpl, 0, aux3, nspin,
889+ & DRho_gga, Ex, Ec, DEx, DEc, Vscf_gga,
890+ & dummy_DVxcdn, stressl )
891+
892+ call timer("BSC-CellXC",2)
893+
894+ else
895+
896+ call timer("GXC-CellXC",1)
897+
898+ ! Note that RG's meshLim is in blocks of nsm points, and 1-based
899+ ! Example.
900+ ! 1:16 17:32 (base 1)
901+ ! 0:15 16:31 (base 0)
902+ ! Times nsm=2:
903+ ! 0:30 32:62
904+ ! Add 1 to lb, and nsm to ub:
905+ ! 1:32 33:64 (base 1, in units of small-points)
906+
907+ myBox(1,:) = (meshLim(1,:)-1)*nsm + 1
908+ myBox(2,:) = (meshLim(2,:)-1)*nsm + nsm
909+
910+ call cellXC( 0, cell, ntm, myBox(1,1), myBox(2,1),
911+ & myBox(1,2), myBox(2,2),
912+ & myBox(1,3), myBox(2,3), nspin,
913+ & DRho_gga, Ex, Ec, DEx, DEc, stressXC, Vscf_gga,
914+ & keep_input_distribution = .true. )
915+ ! Vscf is still sequential after the call to JMS's cellxc
916+ stress = stress + stressXC
917+ call timer("GXC-CellXC",2)
918+
919+ endif ! use_bsc_cellxc
920+
921+ ! Go back to uniform distribution
922+ if (nodes.gt.1) then
923+ call setMeshDistr( UNIFORM, nsm, nsp, nml, nmpl, ntml, ntpl)
924+ endif
925+
926+ ! Redistribute to the Vxc array
927+ do ispin = 1,nspin
928+ fsrc => Vscf_gga(:,ispin)
929+ fdst => Vscf(:,ispin)
930+ call distMeshData( LINEAR, fsrc, UNIFORM, fdst, KEEP )
931+ enddo
932+ call de_alloc( DRho_gga, 'DRho_gga', 'dhscf' )
933+ call de_alloc( Vscf_gga, 'Vscf_gga', 'dhscf' )
934+
935+ call timer("XC",2)
936
937 if ( debug_dhscf ) then
938 write(*,debug_fmt) Node,'XC',
939 & (sqrt(sum(Vscf(:,ispin)**2)),ispin=1,nspin)
940 end if
941
942-
943-#ifndef BSC_CELLXC
944-! Vscf is still sequential after the call to JMS's cellxc
945-#else /* BSC_CELLXC */
946- call de_alloc( DRho_gga, 'DRho_gga', 'dhscf' )
947- call de_alloc( Vscf_gga, 'Vscf_gga', 'dhscf' )
948-#endif /* BSC_CELLXC */
949-
950 Exc = Ex + Ec
951 Dxc = DEx + DEc
952
953- call timer("CellXC",2)
954-
955 ! Vscf contains only Vxc, and is UNIFORM and sequential
956 ! Now we add up the other contributions to it, at
957 ! the same time that we get DRho back to true DeltaRho form
958@@ -1618,10 +1601,6 @@
959 enddo
960 !$OMP end parallel
961
962-#ifndef BSC_CELLXC
963- stress = stress + stressXC
964-#endif /* ! BSC_CELLXC */
965-
966 C ----------------------------------------------------------------------
967 C Save total potential
968 C ----------------------------------------------------------------------
969@@ -1700,20 +1679,16 @@
970
971 #ifdef MPI
972 C Global reduction of Uscf/DUscf/Uatm/Enaatm/Enascf
973-#ifndef BSC_CELLXC
974-! Note that Exc and Dxc are already reduced in the new cellxc
975-#endif /* ! BSC_CELLXC */
976 sbuffer(1) = Uscf
977 sbuffer(2) = DUscf
978 sbuffer(3) = Uatm
979 sbuffer(4) = Enaatm
980 sbuffer(5) = Enascf
981-#ifdef BSC_CELLXC
982- sbuffer(6) = Exc
983- sbuffer(7) = Dxc
984-#else
985- sbuffer(6:7) = 0._dp
986-#endif /* BSC_CELLXC */
987+ if (use_bsc_cellxc) then
988+ sbuffer(6) = Exc
989+ sbuffer(7) = Dxc
990+ endif
991+
992 call MPI_AllReduce( sbuffer, rbuffer, 7, MPI_double_precision,
993 & MPI_Sum, MPI_Comm_World, MPIerror )
994 Uscf = rbuffer(1)
995@@ -1721,10 +1696,10 @@
996 Uatm = rbuffer(3)
997 Enaatm = rbuffer(4)
998 Enascf = rbuffer(5)
999-#ifdef BSC_CELLXC
1000- Exc = rbuffer(6)
1001- Dxc = rbuffer(7)
1002-#endif /* BSC_CELLXC */
1003+ if (use_bsc_cellxc) then
1004+ Exc = rbuffer(6)
1005+ Dxc = rbuffer(7)
1006+ endif
1007 #endif /* MPI */
1008
1009 C Add contribution to stress from the derivative of the Jacobian of ---
1010@@ -1782,11 +1757,8 @@
1011 if ( harrisfun) then
1012 ! Forhar deals internally with its own needs
1013 ! for distribution changes
1014-#ifndef BSC_CELLXC
1015+! Upon entry, everything is UNIFORM, sequential form
1016 call forhar( ntpl, nspin, nml, ntml, ntm, npcc, cell,
1017-#else /* BSC_CELLXC */
1018- call forhar( ntpl, nspin, nml, ntml, npcc, cell,
1019-#endif /* BSC_CELLXC */
1020 & rhoatm, rhopcc, Vna, DRho, Vscf, Vaux )
1021 ! Upon return, everything is UNIFORM, sequential form
1022 endif
1023@@ -1993,4 +1965,30 @@
1024 endif
1025 end subroutine delk_wrapper
1026
1027+ !> Chech whether we need to initialize stencils for
1028+ !> BSC's version of cellxc, and whether we have a vdw
1029+ !> functional
1030+ subroutine xc_setup(need_grads,is_vdw)
1031+ use siestaxc, only: getxc
1032+
1033+ logical, intent(out) :: need_grads
1034+ logical, intent(out) :: is_vdw
1035+
1036+ integer :: nf, nXCfunc
1037+ character(len=20):: XCauth(10), XCfunc(10)
1038+
1039+ need_grads = .false.
1040+ is_vdw = .false.
1041+ call getXC( nXCfunc, XCfunc, XCauth )
1042+ do nf = 1,nXCfunc
1043+ if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then
1044+ need_grads = .true.
1045+ endif
1046+ if ( XCfunc(nf).eq.'VDW' .or. XCfunc(nf).eq.'vdw') then
1047+ need_grads = .true.
1048+ is_vdw = .true.
1049+ endif
1050+ enddo
1051+ end subroutine xc_setup
1052+
1053 end module m_dhscf
1054
1055=== modified file 'Src/forhar.F'
1056--- Src/forhar.F 2016-01-25 16:00:16 +0000
1057+++ Src/forhar.F 2019-01-30 13:43:54 +0000
1058@@ -6,35 +6,15 @@
1059 ! See Docs/Contributors.txt for a list of contributors.
1060 !
1061 module m_forhar
1062- use precision, only : dp, grid_p
1063- use alloc, only : re_alloc, de_alloc
1064-#ifndef BSC_CELLXC
1065- use parallel, only : ProcessorY
1066- use mesh, only : NSM
1067- use siestaXC, only : cellXC ! Finds xc energy and potential
1068- use siestaXC, only : myMeshBox ! Returns my processor mesh box
1069- use siestaXC, only : jms_setMeshDistr => setMeshDistr
1070- ! Sets a distribution of mesh
1071- ! points over parallel processors
1072
1073-#else /* BSC_CELLXC */
1074- use mesh
1075- use moreMeshSubs, only : setMeshDistr, distMeshData
1076- use moreMeshSubs, only: UNIFORM, LINEAR
1077- use moreMeshSubs, only: KEEP
1078-#endif /* BSC_CELLXC */
1079 implicit none
1080
1081 public :: forhar
1082 private
1083
1084 CONTAINS
1085-#ifndef BSC_CELLXC
1086 subroutine forhar( NTPL, NSPIN, NML, NTML, NTM, NPCC,
1087 $ CELL, RHOATM,
1088-#else /* BSC_CELLXC */
1089- subroutine forhar( NTPL, NSPIN, NML, NTML, NPCC, CELL, RHOATM,
1090-#endif /* BSC_CELLXC */
1091 & RHOPCC, VNA, DRHOOUT, VHARRIS1, VHARRIS2 )
1092
1093 C **********************************************************************
1094@@ -45,6 +25,7 @@
1095 C In the first SCF step, V_Hartree(DeltaRho_in) is zero, because
1096 C in that case, Rho_SCF(r) = Rho_atm(r) and therefore, DeltaRho(r) = 0
1097 C This calculation will be skipped.
1098+C NOTE: This is true only if the initial DM is built from atomic charges...
1099 C If Harris + Spin polarized in the first SCF step, then Vharris2 will
1100 C multiply to D Rho(Harris)/D R inside dfscf, and the change of
1101 C the harris density respect the displacement
1102@@ -53,25 +34,30 @@
1103 C Coded by J. Junquera 09/00
1104 C **********************************************************************
1105
1106+ use precision, only : dp, grid_p
1107+ use alloc, only : re_alloc, de_alloc
1108+
1109+ use parallel, only : nodes
1110+ use mesh, only : NSM, nsp, meshLim
1111+ use siestaXC, only : cellXC ! Finds xc energy and potential
1112+
1113+ use moreMeshSubs, only : setMeshDistr, distMeshData
1114+ use moreMeshSubs, only: UNIFORM, LINEAR
1115+ use moreMeshSubs, only: KEEP
1116+
1117+ use fdf, only: fdf_get
1118+
1119 INTEGER :: NTPL, NML(3), NTML(3)
1120-#ifndef BSC_CELLXC
1121 INTEGER, INTENT(IN) :: NSPIN, NPCC, NTM(3)
1122-#else /* BSC_CELLXC */
1123- INTEGER, INTENT(IN) :: NSPIN, NPCC
1124-#endif /* BSC_CELLXC */
1125
1126 REAL(dp), INTENT(IN) :: CELL(3,3)
1127 REAL(grid_p), INTENT(IN) :: VNA(NTPL), RHOATM(NTPL),
1128 & RHOPCC(NTPL)
1129- REAL(grid_p), INTENT(INOUT) :: DRHOOUT(NTPL,NSPIN)
1130- REAL(grid_p), TARGET, INTENT(INOUT) :: VHARRIS1(NTPL,NSPIN)
1131- REAL(grid_p), INTENT(INOUT) :: VHARRIS2(NTPL)
1132+ REAL(grid_p), INTENT(IN) :: DRHOOUT(NTPL,NSPIN)
1133+ REAL(grid_p), TARGET, INTENT(OUT) :: VHARRIS1(NTPL,NSPIN)
1134+ REAL(grid_p), INTENT(OUT) :: VHARRIS2(NTPL)
1135
1136-#ifndef BSC_CELLXC
1137- EXTERNAL REORD
1138-#else /* BSC_CELLXC */
1139 EXTERNAL bsc_cellxc
1140-#endif /* BSC_CELLXC */
1141
1142 ! AG: Note: REAL*4 variables are really REAL(kind=grid_p)
1143 !
1144@@ -92,7 +78,7 @@
1145 C is Drho_out-Rhoatm.
1146 C ***** OUTPUT *********************************************************
1147 C REAL*4 VHARRIS1(NTPL,NSPIN) : Vna + V_Hartree(DeltaRho_in) + V_xc(Rho_in)
1148-C REAL*4 VHARRIS2(NTPL) : Vna + V_Hartree(Rho_in) +
1149+C REAL*4 VHARRIS2(NTPL) : Vna + V_Hartree(DeltaRho_in) +
1150 C - DV_xc(Rho_in)/DRho_in * (Rho_out-Rho_in)
1151 C If Harris forces are computed in the
1152 C first SCF step, it does not depend on spin
1153@@ -104,24 +90,16 @@
1154 C ----------------------------------------------------------------------
1155 C Internal variables and arrays
1156 C ----------------------------------------------------------------------
1157-#ifndef BSC_CELLXC
1158- INTEGER IP, ISPIN, ISPIN2, myBox(2,3)
1159+
1160+ INTEGER IP, ISPIN, ISPIN2, myBox(2,3), NMPL
1161 REAL(dp) EX, EC, DEX, DEC, STRESS(3,3)
1162- INTEGER, SAVE:: JDGdistr=-1
1163
1164- real(grid_p), pointer :: drhoin(:,:),
1165- & dvxcdn(:,:,:)
1166-#else /* BSC_CELLXC */
1167- INTEGER :: IP, ISPIN, ISPIN2, NMPL
1168- REAL(dp) :: EX, EC, DEX, DEC, STRESSL(3,3)
1169 real(grid_p) :: aux3(3,1) !! dummy arrays for cellxc
1170-
1171 real(grid_p), pointer :: drhoin(:,:), drhoin_par(:,:),
1172 & dvxcdn(:,:,:), dvxcdn_par(:,:,:),
1173 & vharris1_par(:,:), fsrc(:), fdst(:)
1174- INTEGER :: ntpl_3
1175-#endif /* BSC_CELLXC */
1176-
1177+ logical :: use_bsc_cellxc
1178+
1179 nullify( drhoin, dvxcdn )
1180 call re_alloc( drhoin, 1, ntpl, 1, nspin, 'drhoin', 'forhar' )
1181 call re_alloc( dvxcdn, 1, ntpl, 1, nspin, 1, nspin,
1182@@ -134,18 +112,6 @@
1183 VHARRIS2(:) = 0.0_grid_p
1184 DRHOIN(:,:) = 0.0_grid_p
1185 DVXCDN(:,:,:) = 0.0_grid_p
1186-#ifndef BSC_CELLXC
1187-
1188-C ----------------------------------------------------------------------
1189-C Set uniform distribution of mesh points and find my processor mesh box
1190-C ----------------------------------------------------------------------
1191-
1192- call jms_setMeshDistr( distrID=JDGdistr, nMesh=NTM,
1193- . nNodesX=1, nNodesY=ProcessorY, nBlock=NSM )
1194- call myMeshBox( NTM, JDGdistr, myBox )
1195-#else /* BSC_CELLXC */
1196- STRESSL(:,:) = 0.0_dp
1197-#endif /* BSC_CELLXC */
1198
1199 C ----------------------------------------------------------------------
1200 C Compute exchange-correlation energy and potential and
1201@@ -153,6 +119,7 @@
1202 C or the sum of atomic charges.
1203 C ----------------------------------------------------------------------
1204
1205+ ! All these arrays are in the UNIFORM distribution,in SEQUENTIAL form
1206 DO ISPIN = 1, NSPIN
1207 DRHOIN(1:NTPL,ISPIN) = RHOATM(1:NTPL)/NSPIN
1208 IF (NPCC .EQ. 1)
1209@@ -160,79 +127,71 @@
1210 . RHOPCC(1:NTPL)/NSPIN
1211 ENDDO
1212
1213-#ifdef BSC_CELLXC
1214-#ifdef MPI
1215- ! The input distribution is UNIFORM, but we need to work
1216- ! with the LINEAR one
1217- call setMeshDistr( LINEAR, nsm, nsp, nml, nmpl, ntml, ntpl )
1218-#endif
1219- ntpl_3 = ntpl
1220+ ! Give the opportunity to use BSC's version
1221+ use_bsc_cellxc = fdf_get("XC.Use.BSC.Cellxc",.false.)
1222+
1223+ ! The input distribution is UNIFORM, but we need to work with the
1224+ ! "zero/not-zero rho" distribution (miscalled 'LINEAR')
1225+ if (nodes.gt.1) then
1226+ call setMeshDistr( LINEAR, nsm, nsp, nml, nmpl, ntml, ntpl )
1227+ endif
1228
1229 nullify( drhoin_par, vharris1_par, dvxcdn_par )
1230- call re_alloc( drhoin_par, 1, ntpl_3, 1, nspin,
1231+ call re_alloc( drhoin_par, 1, ntpl, 1, nspin,
1232 & 'drhoin_par', 'forhar' )
1233- call re_alloc( vharris1_par, 1, ntpl_3, 1, nspin,
1234+ call re_alloc( vharris1_par, 1, ntpl, 1, nspin,
1235 & 'vharris1_par', 'forhar' )
1236- call re_alloc( dvxcdn_par, 1, ntpl_3, 1, nspin, 1, nspin,
1237+ call re_alloc( dvxcdn_par, 1, ntpl, 1, nspin, 1, nspin,
1238 & 'dvxcdn_par', 'forhar' )
1239-#endif /* BSC_CELLXC */
1240
1241 DO ISPIN = 1, NSPIN
1242-#ifndef BSC_CELLXC
1243- CALL REORD(DRHOIN(1,ISPIN),DRHOIN(1,ISPIN),NML,NSM,+1)
1244-#else /* BSC_CELLXC */
1245 fsrc => drhoin(:,ispin)
1246 fdst => drhoin_par(:,ispin)
1247 call distMeshData( UNIFORM, fsrc, LINEAR, fdst, KEEP )
1248-#endif /* BSC_CELLXC */
1249 ENDDO
1250
1251-#ifndef BSC_CELLXC
1252- CALL CELLXC( 0, CELL, NTM, myBox(1,1), myBox(2,1),
1253- . myBox(1,2), myBox(2,2),
1254- . myBox(1,3), myBox(2,3), NSPIN, DRHOIN,
1255- . EX, EC, DEX, DEC, STRESS, VHARRIS1, DVXCDN )
1256-#else /* BSC_CELLXC */
1257- CALL bsc_cellxc( 0, 1, CELL, NTML, NTML, NTPL, 0, AUX3, NSPIN,
1258+ if (use_bsc_cellxc) then
1259+
1260+ STRESS(:,:) = 0.0_dp
1261+ CALL bsc_cellxc( 0, 1, CELL, NTML, NTML, NTPL, 0, AUX3, NSPIN,
1262 & DRHOIN_PAR, EX, EC, DEX, DEC, VHARRIS1_PAR,
1263- & DVXCDN_PAR, STRESSL )
1264-#endif /* BSC_CELLXC */
1265+ & DVXCDN_PAR, STRESS )
1266+
1267+ else
1268+
1269+ myBox(1,:) = (meshLim(1,:)-1)*nsm + 1
1270+ myBox(2,:) = (meshLim(2,:)-1)*nsm + nsm
1271+
1272+ CALL CELLXC( 0, CELL, NTM, myBox(1,1), myBox(2,1),
1273+ . myBox(1,2), myBox(2,2),
1274+ . myBox(1,3), myBox(2,3), NSPIN, DRHOIN_par,
1275+ . EX, EC, DEX, DEC, STRESS, VHARRIS1_par, DVXCDN_par,
1276+ & keep_input_distribution = .true. )
1277+
1278+ endif
1279+
1280+ if (nodes.gt.1) then
1281+! Everything back to UNIFORM, sequential
1282+ call setMeshDistr( UNIFORM, nsm, nsp,
1283+ & nml, nmpl, ntml, ntpl )
1284+ endif
1285
1286 DO ISPIN = 1, NSPIN
1287-#ifndef BSC_CELLXC
1288- CALL REORD(DRHOIN(1,ISPIN),DRHOIN(1,ISPIN),NML,NSM,-1)
1289- CALL REORD(VHARRIS1(1,ISPIN),VHARRIS1(1,ISPIN),NML,NSM,-1)
1290-
1291-#else /* BSC_CELLXC */
1292 fsrc => VHARRIS1_PAR(:,ISPIN)
1293 fdst => VHARRIS1(:,ispin)
1294 call distMeshData( LINEAR, fsrc, UNIFORM, fdst, KEEP )
1295-#endif /* BSC_CELLXC */
1296 DO ISPIN2 = 1, NSPIN
1297-#ifndef BSC_CELLXC
1298- CALL REORD(DVXCDN(1,ISPIN,ISPIN2),DVXCDN(1,ISPIN,ISPIN2),
1299- . NML,NSM,-1)
1300-#else /* BSC_CELLXC */
1301 fsrc => DVXCDN_PAR(:,ISPIN,ISPIN2)
1302 fdst => DVXCDN(:,ISPIN,ISPIN2)
1303 call distMeshData( LINEAR, fsrc, UNIFORM, fdst, KEEP )
1304-#endif /* BSC_CELLXC */
1305 ENDDO
1306 ENDDO
1307-
1308-#ifdef BSC_CELLXC
1309 call de_alloc( dvxcdn_par, 'dvxcdn_par', 'forhar' )
1310 call de_alloc( vharris1_par, 'vharris1_par', 'forhar' )
1311 call de_alloc( drhoin_par, 'drhoin_par', 'forhar' )
1312
1313-#ifdef MPI
1314- call setMeshDistr( UNIFORM, nsm, nsp, nml, nmpl, ntml, ntpl )
1315-#endif
1316-#endif /* BSC_CELLXC */
1317+
1318 DO ISPIN = 1, NSPIN
1319- IF( NPCC .EQ. 1)
1320- & DRHOIN(1:NTPL,ISPIN) = DRHOIN(1:NTPL,ISPIN) -
1321- & RHOPCC(1:NTPL)/NSPIN
1322 DO IP = 1, NTPL
1323 VHARRIS1(IP,ISPIN) = VHARRIS1(IP,ISPIN) + VNA(IP)
1324 ENDDO
1325@@ -255,7 +214,7 @@
1326
1327 C ----------------------------------------------------------------------
1328 C Since V_Hartree(DeltaRho_in) = 0.0, we only add to vharris2 the neutral
1329-C atom potential
1330+C atom potential (note sign change of above intermediate result)
1331 C ----------------------------------------------------------------------
1332 DO IP = 1, NTPL
1333 VHARRIS2(IP) = VNA(IP) - VHARRIS2(IP)
1334
1335=== modified file 'Src/ldau.F'
1336--- Src/ldau.F 2016-06-22 19:57:38 +0000
1337+++ Src/ldau.F 2019-01-30 13:43:54 +0000
1338@@ -197,14 +197,6 @@
1339 type(rad_func), pointer :: pp
1340 C ......................
1341
1342-#ifdef BSC_CELLXC
1343- ! The only reason is that the XC functional
1344- ! is not calculated exactly equivalent.
1345- ! As such, one may delete this below line to make them accessible
1346- ! in conjunction with each other...
1347- call die('LDA+U and BSC_CELLXC is not compatible')
1348-#endif
1349-
1350 ! Start time counter
1351 call timer( 'hubbard_term', 1 )
1352
1353
1354=== modified file 'Src/meshsubs.F'
1355--- Src/meshsubs.F 2018-12-19 20:31:23 +0000
1356+++ Src/meshsubs.F 2019-01-30 13:43:54 +0000
1357@@ -1323,7 +1323,8 @@
1358
1359 !------------------------------------------------------------------
1360 subroutine distriPhiOnMesh( nm, nmpl, norb, iaorb,
1361- & iphorb, isa, numphi )
1362+ & iphorb, isa, numphi,
1363+ & need_gradients_in_xc )
1364
1365 C Computes the number of orbitals that intersects every mesh point
1366 C and calls initMeshDistr in order to compute a new data distribution
1367@@ -1337,7 +1338,7 @@
1368 C integer iaorb(norb) : Atom to which each orbital belongs
1369 C integer iphorb(norb) : Orbital index (within atom) of each orbital
1370 C integer isa(na) : Species index of all atoms in supercell
1371-C
1372+C logical need_gradients_in_xc :
1373 C OUTPUT:
1374 C The output values are in the module moreMeshSubs (New limits of the mesh
1375 C and the information needed to transfer data between data
1376@@ -1377,13 +1378,14 @@
1377 use moreMeshSubs, only: initMeshDistr
1378 use moreMeshSubs, only: allocASynBuffer
1379 use moreMeshSubs, only: UNIFORM, QUADRATIC, LINEAR
1380- use cellxc_mod, only: GGA, setGGA
1381 use alloc, only: re_alloc, de_alloc
1382+
1383 implicit none
1384 C Passed arguments
1385 integer, intent(in) :: nm(3), nmpl, norb, iaorb(norb),
1386 & iphorb(norb), isa(*)
1387 integer, intent(out) :: numphi(nmpl)
1388+ logical, intent(in) :: need_gradients_in_xc
1389 C Local variables
1390 integer :: ii, io, ia, iphi, is, iop, ip0, isp
1391 real(dp) :: r2o, dxsp(3,nsp), r2sp(nsp)
1392@@ -1442,7 +1444,7 @@
1393 endif
1394 enddo
1395 call initMeshDistr( UNIFORM, LINEAR, nm, wload )
1396- if (GGA) call initMeshExtencil( LINEAR, nm )
1397+ if (need_gradients_in_xc) call initMeshExtencil( LINEAR, nm )
1398
1399 C Free local memory
1400 call de_alloc( wload, 'wload', 'distriPhiOnMesh' )
1401
1402=== modified file 'Src/siesta_init.F'
1403--- Src/siesta_init.F 2018-10-15 21:10:25 +0000
1404+++ Src/siesta_init.F 2019-01-30 13:43:54 +0000
1405@@ -30,9 +30,6 @@
1406 & initatomlists, no_l, rmaxldau
1407 use fdf
1408 use sys, only: die, bye
1409-#ifdef BSC_CELLXC
1410- use bsc_xcmod, only: setXC
1411-#endif /* BSC_CELLXC */
1412 use molecularmechanics, only: inittwobody
1413 use metaforce, only: initmeta
1414 use m_mpi_utils, only: broadcast
1415@@ -366,11 +363,7 @@
1416
1417
1418 ! Initialise exchange-correlation functional information
1419-#ifndef BSC_CELLXC
1420 call read_xc_info()
1421-#else /* BSC_CELLXC */
1422- call setXC( )
1423-#endif /* BSC_CELLXC */
1424
1425 ! Initialise force field component
1426 call inittwobody( )
1427
1428=== removed file 'Src/xc.f'
1429--- Src/xc.f 2017-07-01 20:58:50 +0000
1430+++ Src/xc.f 1970-01-01 00:00:00 +0000
1431@@ -1,2547 +0,0 @@
1432-!
1433-! Copyright (C) 1996-2016 The SIESTA group
1434-! This file is distributed under the terms of the
1435-! GNU General Public License: see COPYING in the top directory
1436-! or http://www.gnu.org/copyleft/gpl.txt.
1437-! See Docs/Contributors.txt for a list of contributors.
1438-!
1439-! *******************************************************************
1440-! This file contains XC subroutines used when siesta is compiled with
1441-! option BSC_CELLXC. Otherwise, the SiestaXC library is used.
1442-! *******************************************************************
1443-
1444- subroutine atomxc( IREL, NR, MAXR, RMESH, nspin, Dens,
1445- . EX, EC, DX, DC, VXC )
1446-
1447-C *******************************************************************
1448-C Finds total exchange-correlation energy and potential for a
1449-C spherical electron density distribution.
1450-C This version implements the Local (spin) Density Approximation and
1451-C the Generalized-Gradient-Aproximation with the 'explicit mesh
1452-C functional' approach of White & Bird, PRB 50, 4954 (1994).
1453-C Gradients are 'defined' by numerical derivatives, using 2*NN+1 mesh
1454-C points, where NN is a parameter defined below
1455-C Ref: L.C.Balbas et al, PRB 64, 165110 (2001)
1456-C Wrtten by J.M.Soler using algorithms developed by
1457-C L.C.Balbas, J.L.Martins and J.M.Soler, Dec.1996
1458-C ************************* INPUT ***********************************
1459-C CHARACTER*(*) FUNCTL : Functional to be used:
1460-C 'LDA' or 'LSD' => Local (spin) Density Approximation
1461-C 'GGA' => Generalized Gradient Corrections
1462-C Uppercase is optional
1463-C CHARACTER*(*) AUTHOR : Parametrization desired:
1464-C 'CA' or 'PZ' => LSD Perdew & Zunger, PRB 23, 5075 (1981)
1465-C 'PW91' => GGA Perdew & Wang, JCP, 100, 1290 (1994)
1466-C 'PW92' => LSD Perdew & Wang, PRB, 45, 13244 (1992). This is
1467-C the local density limit of the next:
1468-C 'PBE' => GGA Perdew, Burke & Ernzerhof, PRL 77, 3865 (1996)
1469-C 'RPBE' => GGA Hammer, Hansen & Norskov, PRB 59, 7413 (1999)
1470-C 'REVPBE' => GGA Zhang & Yang, PRL 80,890(1998)
1471-C 'LYP' => GGA Becke-Lee-Yang-Parr (see subroutine blypxc)
1472-C 'WC' => GGA Wu-Cohen (see subroutine wcxc)
1473-C 'PBESOL' => GGA Perdew et al, PRL, 100, 136406 (2008)
1474-C Uppercase is optional
1475-C INTEGER IREL : Relativistic exchange? (0=>no, 1=>yes)
1476-C INTEGER NR : Number of radial mesh points
1477-C INTEGER MAXR : Physical first dimension of RMESH, Dens and VXC
1478-C REAL*8 RMESH(MAXR) : Radial mesh points
1479-C INTEGER nspin : nspin=1 => unpolarized; nspin=2 => polarized
1480-C REAL*8 Dens(MAXR,nspin) : Total (nspin=1) or spin (nspin=2) electron
1481-C density at mesh points
1482-C ************************* OUTPUT **********************************
1483-C REAL*8 EX : Total exchange energy
1484-C REAL*8 EC : Total correlation energy
1485-C REAL*8 DX : IntegralOf( rho * (eps_x - v_x) )
1486-C REAL*8 DC : IntegralOf( rho * (eps_c - v_c) )
1487-C REAL*8 VXC(MAXR,nspin) : (Spin) exch-corr potential
1488-C ************************ UNITS ************************************
1489-C Distances in atomic units (Bohr).
1490-C Densities in atomic units (electrons/Bohr**3)
1491-C Energy unit depending of parameter EUNIT below
1492-C ********* ROUTINES CALLED *****************************************
1493-C GGAXC, LDAXC
1494-C *******************************************************************
1495-
1496- use precision, only : dp
1497- use bsc_xcmod, only : nXCfunc, XCfunc, XCauth
1498- use bsc_xcmod, only : XCweightX, XCweightC
1499- use sys, only: die
1500- use alloc, only: re_alloc, de_alloc
1501-
1502-C Next line is nonstandard but may be suppressed
1503- implicit none
1504-
1505-C Argument types and dimensions
1506- integer, intent(in) :: IREL
1507- integer, intent(in) :: MAXR
1508- integer, intent(in) :: NR
1509- integer, intent(in) :: nspin
1510- real(dp), intent(in) :: Dens(MAXR,nspin)
1511- real(dp), intent(in) :: RMESH(MAXR)
1512- real(dp), intent(out) :: VXC(MAXR,nspin)
1513- real(dp), intent(out) :: DC
1514- real(dp), intent(out) :: DX
1515- real(dp), intent(out) :: EC
1516- real(dp), intent(out) :: EX
1517-
1518-C Internal parameters
1519-C NN : order of the numerical derivatives: the number of radial
1520-C points used is 2*NN+1
1521-C mspin : must be equal or larger than nspin (4 for non-collinear spin)
1522- integer, parameter :: mspin = 4
1523- integer, parameter :: NN = 5
1524-
1525-C Fix energy unit: EUNIT=1.0 => Hartrees,
1526-C EUNIT=0.5 => Rydbergs,
1527-C EUNIT=0.03674903 => eV
1528- real(dp), parameter :: EUNIT = 0.5_dp
1529-
1530-C DVMIN is added to differential of volume to avoid division by zero
1531- real(dp), parameter :: DVMIN = 1.0d-12
1532-
1533-C Local variables and arrays
1534- logical
1535- . GGA, GGAfunc
1536- integer
1537- . IN, IN1, IN2, IR, IS, JN, NF
1538- real(dp)
1539- . D(mspin), DECDD(mspin), DECDGD(3,mspin),
1540- . DEXDD(mspin), DEXDGD(3,mspin),
1541- . DGDM(-NN:NN), DGIDFJ(-NN:NN), DRDM, DVol,
1542- . DVCDN(mspin,mspin), DVXDN(mspin,mspin),
1543- . EPSC, EPSX, F1, F2, GD(3,mspin), PI
1544- real(dp), pointer :: Aux(:)
1545- external
1546- . GGAXC, LDAXC
1547-
1548-C Set GGA switch
1549- GGA = .false.
1550- do nf = 1,nXCfunc
1551- if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then
1552- GGA = .true.
1553- else
1554- if ( XCfunc(nf).ne.'LDA' .and. XCfunc(nf).ne.'lda' .and.
1555- . XCfunc(nf).ne.'LSD' .and. XCfunc(nf).ne.'lsd' ) then
1556- call die('ATOMXC: Unknown functional ' // XCfunc(nf))
1557- endif
1558- endif
1559- enddo
1560-
1561-C Initialize output
1562- EX = 0.0_dp
1563- EC = 0.0_dp
1564- DX = 0.0_dp
1565- DC = 0.0_dp
1566- do IS = 1,nspin
1567- do IR = 1,NR
1568- VXC(IR,IS) = 0.0_dp
1569- enddo
1570- enddo
1571-
1572-C Set up workspace array
1573- if (GGA) then
1574- nullify( Aux )
1575- call re_alloc( AUX, 1, NR, 'AUX', 'atomxc' )
1576- endif
1577-
1578-C Get number pi
1579- PI = 4.0_dp * ATAN(1.0_dp)
1580-
1581-C Loop on mesh points
1582- do IR = 1,NR
1583-
1584-C Find interval of neighbour points to calculate derivatives
1585- IN1 = MAX( 1, IR-NN ) - IR
1586- IN2 = MIN( NR, IR+NN ) - IR
1587-
1588-C Find weights of numerical derivation from Lagrange
1589-C interpolation formula
1590- do IN = IN1,IN2
1591- IF (IN .EQ. 0) THEN
1592- DGDM(IN) = 0
1593- do JN = IN1,IN2
1594- IF (JN.NE.0) DGDM(IN) = DGDM(IN) + 1.D0 / (0 - JN)
1595- enddo
1596- ELSE
1597- F1 = 1
1598- F2 = 1
1599- do JN = IN1,IN2
1600- IF (JN.NE.IN .AND. JN.NE.0) F1 = F1 * (0 - JN)
1601- IF (JN.NE.IN) F2 = F2 * (IN - JN)
1602- enddo
1603- DGDM(IN) = F1 / F2
1604- ENDIF
1605- enddo
1606-
1607-C Find dr/dmesh
1608- DRDM = 0.0_dp
1609- do IN = IN1,IN2
1610- DRDM = DRDM + RMESH(IR+IN) * DGDM(IN)
1611- enddo
1612-
1613-C Find differential of volume. Use trapezoidal integration rule
1614- DVol = 4.0_dp * PI * RMESH(IR)**2 * DRDM
1615-C DVMIN is a small number added to avoid a division by zero
1616- DVol = DVol + DVMIN
1617- if (IR.eq.1 .or. IR.eq.NR) DVol = 0.5_dp*DVol
1618- if (GGA) Aux(IR) = DVol
1619-
1620-C Find the weights for the derivative d(gradF(i))/d(F(j)), of
1621-C the gradient at point i with respect to the value at point j
1622- if (GGA) then
1623- do IN = IN1,IN2
1624- DGIDFJ(IN) = DGDM(IN) / DRDM
1625- enddo
1626- endif
1627-
1628-C Find density and gradient of density at this point
1629- do IS = 1,nspin
1630- D(IS) = Dens(IR,IS)
1631- enddo
1632- if (GGA) then
1633- do IS = 1,nspin
1634- GD(1,IS) = 0.0_dp
1635- GD(2,IS) = 0.0_dp
1636- GD(3,IS) = 0.0_dp
1637- do IN = IN1,IN2
1638- GD(3,IS) = GD(3,IS) + DGIDFJ(IN) * Dens(IR+IN,IS)
1639- enddo
1640- enddo
1641- endif
1642-
1643-C Loop over exchange-correlation functions
1644- do nf = 1,nXCfunc
1645-
1646-C Is this a GGA?
1647- if (XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then
1648- GGAfunc = .true.
1649- else
1650- GGAfunc = .false.
1651- endif
1652-
1653-C Find exchange and correlation energy densities and their
1654-C derivatives with respect to density and density gradient
1655- if (GGAfunc) then
1656- CALL GGAXC( XCauth(nf), IREL, nspin, D, GD,
1657- . EPSX, EPSC, DEXDD, DECDD, DEXDGD, DECDGD )
1658- else
1659- CALL LDAXC( XCauth(nf), IREL, nspin, D, EPSX, EPSC, DEXDD,
1660- . DECDD, DVXDN, DVCDN )
1661- endif
1662-
1663-C Scale terms by weights
1664- EPSX = EPSX*XCweightX(nf)
1665- EPSC = EPSC*XCweightC(nf)
1666- DEXDD(1:nspin) = DEXDD(1:nspin)*XCweightX(nf)
1667- DECDD(1:nspin) = DECDD(1:nspin)*XCweightC(nf)
1668- if (GGAfunc) then
1669- DEXDGD(1:3,1:nspin) = DEXDGD(1:3,1:nspin)*XCweightX(nf)
1670- DECDGD(1:3,1:nspin) = DECDGD(1:3,1:nspin)*XCweightC(nf)
1671- endif
1672-
1673-C Add contributions to exchange-correlation energy and its
1674-C derivatives with respect to density at all points
1675- do IS = 1,nspin
1676- EX = EX + DVol*D(IS)*EPSX
1677- EC = EC + DVol*D(IS)*EPSC
1678- DX = DX + DVol*D(IS)*(EPSX - DEXDD(IS))
1679- DC = DC + DVol*D(IS)*(EPSC - DECDD(IS))
1680- if (GGAfunc) then
1681- VXC(IR,IS) = VXC(IR,IS) + DVol*(DEXDD(IS) + DECDD(IS))
1682- do IN = IN1,IN2
1683- DX= DX - DVol*Dens(IR+IN,IS)*DEXDGD(3,IS)*DGIDFJ(IN)
1684- DC= DC - DVol*Dens(IR+IN,IS)*DECDGD(3,IS)*DGIDFJ(IN)
1685- VXC(IR+IN,IS) = VXC(IR+IN,IS) +
1686- . DVol*(DEXDGD(3,IS) + DECDGD(3,IS))*DGIDFJ(IN)
1687- enddo
1688- else
1689- if (GGA) then
1690- VXC(IR,IS) = VXC(IR,IS) + DVol*(DEXDD(IS) + DECDD(IS))
1691- else
1692- VXC(IR,IS) = VXC(IR,IS) + DEXDD(IS) + DECDD(IS)
1693- endif
1694- endif
1695- enddo
1696-
1697- enddo
1698-
1699- enddo
1700-
1701-C Divide by volume element to obtain the potential (per electron)
1702- if (GGA) then
1703- do IS = 1,NSPIN
1704- do IR = 1,NR
1705- DVol = AUX(IR)
1706- VXC(IR,IS) = VXC(IR,IS) / DVol
1707- enddo
1708- enddo
1709- call de_alloc( AUX, 'AUX', 'atomxc' )
1710- endif
1711-
1712-C Divide by energy unit
1713- EX = EX / EUNIT
1714- EC = EC / EUNIT
1715- DX = DX / EUNIT
1716- DC = DC / EUNIT
1717- do IS = 1,nspin
1718- do IR = 1,NR
1719- VXC(IR,IS) = VXC(IR,IS) / EUNIT
1720- enddo
1721- enddo
1722-
1723- end
1724-
1725-
1726- subroutine exchng( IREL, NSP, DS, EX, VX )
1727-
1728-C *****************************************************************
1729-C Finds local exchange energy density and potential
1730-C Adapted by J.M.Soler from routine velect of Froyen's
1731-C pseudopotential generation program. Madrid, Jan'97. Version 0.5.
1732-C **** Input ******************************************************
1733-C INTEGER IREL : relativistic-exchange switch (0=no, 1=yes)
1734-C INTEGER NSP : spin-polarizations (1=>unpolarized, 2=>polarized)
1735-C REAL*8 DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
1736-C **** Output *****************************************************
1737-C REAL*8 EX : exchange energy density
1738-C REAL*8 VX(NSP) : (spin-dependent) exchange potential
1739-C **** Units ******************************************************
1740-C Densities in electrons/Bohr**3
1741-C Energies in Hartrees
1742-C *****************************************************************
1743-
1744- use precision, only: dp
1745- implicit none
1746-
1747- integer, intent(in) :: nsp, irel
1748- real(dp), intent(in) :: DS(NSP)
1749- real(dp), intent(out) :: VX(NSP)
1750- real(dp), intent(out) :: EX
1751-
1752- real(dp), parameter :: zero = 0.0_dp, one = 1.0_dp
1753- real(dp), parameter :: pfive = 0.5_dp, opf = 1.5_dp
1754- real(dp), parameter :: C014 = 0.014_dp
1755-
1756- real(dp) :: pi, trd, ftrd, tftm, a0
1757- real(dp) :: alp, d1, d2, d, z, fz, fzp, rs, vxp, exp_var
1758- real(dp) :: beta, sb, vxf, exf, alb
1759-
1760- PI=4*ATAN(ONE)
1761- TRD = ONE/3
1762- FTRD = 4*TRD
1763- TFTM = 2**FTRD-2
1764- A0 = (4/(9*PI))**TRD
1765-
1766-C X-alpha parameter:
1767- ALP = 2 * TRD
1768-
1769- IF (NSP .EQ. 2) THEN
1770- D1 = MAX(DS(1),0.D0)
1771- D2 = MAX(DS(2),0.D0)
1772- D = D1 + D2
1773- IF (D .LE. ZERO) THEN
1774- EX = ZERO
1775- VX(1) = ZERO
1776- VX(2) = ZERO
1777- RETURN
1778- ENDIF
1779- Z = (D1 - D2) / D
1780- FZ = ((1+Z)**FTRD+(1-Z)**FTRD-2)/TFTM
1781- FZP = FTRD*((1+Z)**TRD-(1-Z)**TRD)/TFTM
1782- ELSE
1783- D = DS(1)
1784- IF (D .LE. ZERO) THEN
1785- EX = ZERO
1786- VX(1) = ZERO
1787- RETURN
1788- ENDIF
1789- Z = ZERO
1790- FZ = ZERO
1791- FZP = ZERO
1792- ENDIF
1793- RS = (3 / (4*PI*D) )**TRD
1794- VXP = -(3*ALP/(2*PI*A0*RS))
1795- EXP_VAR = 3*VXP/4
1796- IF (IREL .EQ. 1) THEN
1797- BETA = C014/RS
1798- SB = SQRT(1+BETA*BETA)
1799- ALB = LOG(BETA+SB)
1800- VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
1801- EXP_VAR = EXP_VAR * (ONE-OPF*((BETA*SB-ALB)/BETA**2)**2)
1802- ENDIF
1803- VXF = 2**TRD*VXP
1804- EXF = 2**TRD*EXP_VAR
1805- IF (NSP .EQ. 2) THEN
1806- VX(1) = VXP + FZ*(VXF-VXP) + (1-Z)*FZP*(EXF-EXP_VAR)
1807- VX(2) = VXP + FZ*(VXF-VXP) - (1+Z)*FZP*(EXF-EXP_VAR)
1808- EX = EXP_VAR + FZ*(EXF-EXP_VAR)
1809- ELSE
1810- VX(1) = VXP
1811- EX = EXP_VAR
1812- ENDIF
1813- END
1814-
1815- SUBROUTINE GGAXC( AUTHOR, IREL, nspin, D, GD,
1816- . EPSX, EPSC, DEXDD, DECDD, DEXDGD, DECDGD )
1817-
1818-C Finds the exchange and correlation energies at a point, and their
1819-C derivatives with respect to density and density gradient, in the
1820-C Generalized Gradient Correction approximation.
1821-C Lengths in Bohr, energies in Hartrees
1822-C Written by L.C.Balbas and J.M.Soler, Dec'96. Version 0.5.
1823-C Modified by V.M.Garcia-Suarez to include non-collinear spin. June 2002
1824-
1825- use precision, only : dp
1826- use sys, only : die
1827-
1828- implicit none
1829-
1830- CHARACTER*(*) AUTHOR
1831- INTEGER IREL, nspin, NS, IS, IX
1832- real(dp) THETA, PHI, D(nspin), DECDD(nspin),
1833- . DECDGD(3,nspin), DEXDD(nspin), DEXDGD(3,nspin),
1834- . EPSC, EPSX, GD(3,nspin),
1835- . DD(2), DTOT, DPOL,
1836- . GDD(3,2), TINY, DECDN(2), DEXDN(2),
1837- . VPOL, DECDGN(3,2), DEXDGN(3,2),
1838- . C2, S2, ST, CT, CP, SP, dpolz, dpolxy
1839-
1840-
1841- PARAMETER ( TINY = 1.D-12 )
1842-
1843- IF (nspin .EQ. 4) THEN
1844-C Find eigenvalues of density matrix (up and down densities
1845-C along the spin direction)
1846-C Note: D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
1847- NS = 2
1848- DTOT = D(1) + D(2)
1849-
1850-! Explicit calculation of the rotation-matrix elements from
1851-! the entries of D
1852-
1853- dpolz= D(1)-D(2)
1854- dpolxy= 2.0d0*sqrt(D(3)**2+D(4)**2)
1855- DPOL = sqrt( dpolz**2 + dpolxy**2 )
1856- if ( DPOL.gt.1.0d-12 ) then
1857- THETA = atan2(dpolxy,dpolz)
1858- else
1859- THETA = 0.0_dp
1860- endif
1861- C2 = COS(THETA/2)
1862- S2 = SIN(THETA/2)
1863- ST = SIN(THETA)
1864- CT = COS(THETA)
1865- PHI = ATAN2(-D(4),D(3))
1866- CP = COS(PHI)
1867- SP = SIN(PHI)
1868-
1869- DD(1) = 0.5D0 * ( DTOT + DPOL )
1870- DD(2) = 0.5D0 * ( DTOT - DPOL )
1871-
1872-C Find diagonal elements of the gradient
1873- DO IX = 1,3
1874- GDD(IX,1) = GD(IX,1)*C2**2 + GD(IX,2)*S2**2 +
1875- . 2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
1876- GDD(IX,2) = GD(IX,1)*S2**2 + GD(IX,2)*C2**2 -
1877- . 2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP)
1878- ENDDO
1879- ELSE
1880- NS = nspin
1881- DO 20 IS = 1,nspin
1882-cag Avoid negative densities
1883- DD(IS) = max(D(IS),0.0d0)
1884- DO 30 IX = 1,3
1885- GDD(IX,IS) = GD(IX,IS)
1886- 30 CONTINUE
1887- 20 CONTINUE
1888- ENDIF
1889-
1890- IF (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN
1891- CALL PBEXC( IREL, NS, DD, GDD,
1892- . EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
1893-cmvfs
1894- ELSE IF (AUTHOR.EQ.'RPBE'.OR.AUTHOR.EQ.'rpbe') THEN
1895- CALL RPBEXC( IREL, NS, DD, GDD,
1896- . EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
1897-cmvfs
1898- ELSE IF (AUTHOR.EQ.'WC'.OR.AUTHOR.EQ.'wc') THEN
1899- CALL WCXC( IREL, NS, DD, GDD,
1900- . EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
1901-cea
1902- ELSE IF (AUTHOR.EQ.'REVPBE'.OR.AUTHOR.EQ.'revpbe'
1903- . .OR.AUTHOR.EQ.'revPBE') THEN
1904- CALL REVPBEXC( IREL, NS, DD, GDD,
1905- . EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
1906-cag
1907- ELSE IF (AUTHOR.EQ.'LYP'.OR.AUTHOR.EQ.'lyp') THEN
1908- CALL BLYPXC( NS, DD, GDD, EPSX, EPSC, dEXdn, dECdn,
1909- . DEXDGN, DECDGN)
1910-cag
1911- ELSEIF (AUTHOR.EQ.'PW91' .OR. AUTHOR.EQ.'pw91') THEN
1912- CALL PW91XC( IREL, NS, DD, GDD,
1913- . EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
1914-cjdg
1915- ELSEIF (AUTHOR.EQ.'PBESOL'.OR.AUTHOR.EQ.'pbesol'
1916- . .OR.AUTHOR.EQ.'PBEsol') THEN
1917- CALL PBESOLXC( IREL, NS, DD, GDD,
1918- . EPSX, EPSC, DEXDN, DECDN, DEXDGN, DECDGN )
1919- ELSE
1920- call die('GGAXC: Unknown author ' // trim(AUTHOR))
1921- ENDIF
1922-
1923- IF (nspin .EQ. 4) THEN
1924-C Find dE/dD(ispin) = dE/dDup * dDup/dD(ispin) +
1925-C dE/dDdown * dDown/dD(ispin)
1926-C Note convention:
1927-C DEDD(1)=dE/dD11, DEDD(2)=dE/dD22,
1928-C DEDD(3)=Re(dE/dD12)=Re(dE/dD21),
1929-C DEDD(4)=Im(dE/dD12)=-Im(dE/D21)
1930-C
1931- VPOL = (DEXDN(1)-DEXDN(2)) * CT
1932- DEXDD(1) = 0.5D0 * ( DEXDN(1) + DEXDN(2) + VPOL )
1933- DEXDD(2) = 0.5D0 * ( DEXDN(1) + DEXDN(2) - VPOL )
1934- DEXDD(3) = 0.5d0 * (DEXDN(1)-DEXDN(2)) * ST * CP
1935- DEXDD(4) =-0.5d0 * (DEXDN(1)-DEXDN(2)) * ST * SP
1936- VPOL = (DECDN(1)-DECDN(2)) * CT
1937- DECDD(1) = 0.5D0 * ( DECDN(1) + DECDN(2) + VPOL )
1938- DECDD(2) = 0.5D0 * ( DECDN(1) + DECDN(2) - VPOL )
1939- DECDD(3) = 0.5d0 * (DECDN(1)-DECDN(2)) * ST * CP
1940- DECDD(4) =-0.5d0 * (DECDN(1)-DECDN(2)) * ST * SP
1941-C Gradient terms
1942- DO 40 IX = 1,3
1943- DEXDGD(IX,1) = DEXDGN(IX,1)*C2**2 + DEXDGN(IX,2)*S2**2
1944- DEXDGD(IX,2) = DEXDGN(IX,1)*S2**2 + DEXDGN(IX,2)*C2**2
1945- DEXDGD(IX,3) = 0.5D0*(DEXDGN(IX,1) - DEXDGN(IX,2))*ST*CP
1946- DEXDGD(IX,4) =-0.5D0*(DEXDGN(IX,1) - DEXDGN(IX,2))*ST*SP
1947- DECDGD(IX,1) = DECDGN(IX,1)*C2**2 + DECDGN(IX,2)*S2**2
1948- DECDGD(IX,2) = DECDGN(IX,1)*S2**2 + DECDGN(IX,2)*C2**2
1949- DECDGD(IX,3) = 0.5D0*(DECDGN(IX,1) - DECDGN(IX,2))*ST*CP
1950- DECDGD(IX,4) =-0.5D0*(DECDGN(IX,1) - DECDGN(IX,2))*ST*SP
1951- 40 CONTINUE
1952- ELSE
1953- DO 60 IS = 1,nspin
1954- DEXDD(IS) = DEXDN(IS)
1955- DECDD(IS) = DECDN(IS)
1956- DO 50 IX = 1,3
1957- DEXDGD(IX,IS) = DEXDGN(IX,IS)
1958- DECDGD(IX,IS) = DECDGN(IX,IS)
1959- 50 CONTINUE
1960- 60 CONTINUE
1961- ENDIF
1962-
1963- END
1964-
1965-
1966- SUBROUTINE LDAXC( AUTHOR, IREL, nspin, D, EPSX, EPSC, VX, VC,
1967- . DVXDN, DVCDN )
1968-
1969-C ******************************************************************
1970-C Finds the exchange and correlation energies and potentials, in the
1971-C Local (spin) Density Approximation.
1972-C Written by L.C.Balbas and J.M.Soler, Dec'96.
1973-C Non-collinear spin added by J.M.Soler, May'98
1974-C *********** INPUT ************************************************
1975-C CHARACTER*(*) AUTHOR : Parametrization desired:
1976-C 'CA' or 'PZ' => LSD Perdew & Zunger, PRB 23, 5075 (1981)
1977-C 'PW92' => LSD Perdew & Wang, PRB, 45, 13244 (1992)
1978-C Uppercase is optional
1979-C INTEGER IREL : Relativistic exchange? (0=>no, 1=>yes)
1980-C INTEGER nspin : nspin=1 => unpolarized; nspin=2 => polarized;
1981-C nspin=4 => non-collinear polarization
1982-C REAL*8 D(nspin) : Local (spin) density. For non-collinear
1983-C polarization, the density matrix is given by:
1984-C D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
1985-C *********** OUTPUT ***********************************************
1986-C REAL*8 EPSX, EPSC : Exchange and correlation energy densities
1987-C REAL*8 VX(nspin), VC(nspin) : Exchange and correlation potentials,
1988-C defined as dExc/dD(ispin)
1989-C REAL*8 DVXDN(nspin,nspin) : Derivative of exchange potential with
1990-C respect the charge density, defined
1991-C as DVx(spin1)/Dn(spin2)
1992-C REAL*8 DVCDN(nspin,nspin) : Derivative of correlation potential
1993-C respect the charge density, defined
1994-C as DVc(spin1)/Dn(spin2)
1995-C *********** UNITS ************************************************
1996-C Lengths in Bohr, energies in Hartrees
1997-C ******************************************************************
1998-
1999- use precision, only : dp
2000- use sys, only : die
2001-
2002- implicit none
2003-
2004- CHARACTER*(*) AUTHOR
2005- INTEGER IREL, nspin
2006- real(dp) D(nspin), EPSC, EPSX, VX(nspin), VC(nspin),
2007- . DVXDN(nspin,nspin), DVCDN(nspin,nspin)
2008-
2009- INTEGER IS, NS, ISPIN1, ISPIN2
2010- real(dp) DD(2), DPOL, DTOT, TINY, VCD(2), VPOL, VXD(2)
2011-
2012- PARAMETER ( TINY = 1.D-12 )
2013-
2014- IF (nspin .EQ. 4) THEN
2015-C Find eigenvalues of density matrix (up and down densities
2016-C along the spin direction)
2017-C Note: D(1)=D11, D(2)=D22, D(3)=Real(D12), D(4)=Im(D12)
2018- NS = 2
2019- DTOT = D(1) + D(2)
2020- DPOL = SQRT( (D(1)-D(2))**2 + 4.D0*(D(3)**2+D(4)**2) )
2021- DD(1) = 0.5D0 * ( DTOT + DPOL )
2022- DD(2) = 0.5D0 * ( DTOT - DPOL )
2023- ELSE
2024- NS = nspin
2025- DO 10 IS = 1,nspin
2026-cag Avoid negative densities
2027- DD(IS) = max(D(IS),0.0d0)
2028- 10 CONTINUE
2029- ENDIF
2030-
2031-
2032- DO ISPIN2 = 1, nspin
2033- DO ISPIN1 = 1, nspin
2034- DVXDN(ISPIN1,ISPIN2) = 0.D0
2035- DVCDN(ISPIN1,ISPIN2) = 0.D0
2036- ENDDO
2037- ENDDO
2038-
2039- IF ( AUTHOR.EQ.'CA' .OR. AUTHOR.EQ.'ca' .OR.
2040- . AUTHOR.EQ.'PZ' .OR. AUTHOR.EQ.'pz') THEN
2041- CALL PZXC( IREL, NS, DD, EPSX, EPSC, VXD, VCD, DVXDN, DVCDN )
2042- ELSEIF ( AUTHOR.EQ.'PW92' .OR. AUTHOR.EQ.'pw92' ) THEN
2043- CALL PW92XC( IREL, NS, DD, EPSX, EPSC, VXD, VCD )
2044- ELSE
2045- call die('LDAXC: Unknown author ' // trim(AUTHOR))
2046- ENDIF
2047-
2048- IF (nspin .EQ. 4) THEN
2049-C Find dE/dD(ispin) = dE/dDup * dDup/dD(ispin) +
2050-C dE/dDdown * dDown/dD(ispin)
2051- VPOL = (VXD(1)-VXD(2)) * (D(1)-D(2)) / (DPOL+TINY)
2052- VX(1) = 0.5D0 * ( VXD(1) + VXD(2) + VPOL )
2053- VX(2) = 0.5D0 * ( VXD(1) + VXD(2) - VPOL )
2054- VX(3) = (VXD(1)-VXD(2)) * D(3) / (DPOL+TINY)
2055- VX(4) = (VXD(1)-VXD(2)) * D(4) / (DPOL+TINY)
2056- VPOL = (VCD(1)-VCD(2)) * (D(1)-D(2)) / (DPOL+TINY)
2057- VC(1) = 0.5D0 * ( VCD(1) + VCD(2) + VPOL )
2058- VC(2) = 0.5D0 * ( VCD(1) + VCD(2) - VPOL )
2059- VC(3) = (VCD(1)-VCD(2)) * D(3) / (DPOL+TINY)
2060- VC(4) = (VCD(1)-VCD(2)) * D(4) / (DPOL+TINY)
2061- ELSE
2062- DO 20 IS = 1,nspin
2063- VX(IS) = VXD(IS)
2064- VC(IS) = VCD(IS)
2065- 20 CONTINUE
2066- ENDIF
2067- END
2068-
2069-
2070- SUBROUTINE PBEXC( IREL, nspin, Dens, GDens,
2071- . EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
2072-
2073-C *********************************************************************
2074-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
2075-C Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
2076-C Written by L.C.Balbas and J.M.Soler. December 1996. Version 0.5.
2077-C ******** INPUT ******************************************************
2078-C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
2079-C INTEGER nspin : Number of spin polarizations (1 or 2)
2080-C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
2081-C spin electron density (if nspin=2)
2082-C REAL*8 GDens(3,nspin) : Total or spin density gradient
2083-C ******** OUTPUT *****************************************************
2084-C REAL*8 EX : Exchange energy density
2085-C REAL*8 EC : Correlation energy density
2086-C REAL*8 DEXDD(nspin) : Partial derivative
2087-C d(DensTot*Ex)/dDens(ispin),
2088-C where DensTot = Sum_ispin( Dens(ispin) )
2089-C For a constant density, this is the
2090-C exchange potential
2091-C REAL*8 DECDD(nspin) : Partial derivative
2092-C d(DensTot*Ec)/dDens(ispin),
2093-C where DensTot = Sum_ispin( Dens(ispin) )
2094-C For a constant density, this is the
2095-C correlation potential
2096-C REAL*8 DEXDGD(3,nspin): Partial derivative
2097-C d(DensTot*Ex)/d(GradDens(i,ispin))
2098-C REAL*8 DECDGD(3,nspin): Partial derivative
2099-C d(DensTot*Ec)/d(GradDens(i,ispin))
2100-C ********* UNITS ****************************************************
2101-C Lengths in Bohr
2102-C Densities in electrons per Bohr**3
2103-C Energies in Hartrees
2104-C Gradient vectors in cartesian coordinates
2105-C ********* ROUTINES CALLED ******************************************
2106-C EXCHNG, PW92C
2107-C ********************************************************************
2108-
2109- use precision, only : dp
2110-
2111- implicit none
2112- INTEGER IREL, nspin
2113- real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
2114- . DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
2115-
2116-C Internal variables
2117- INTEGER
2118- . IS, IX
2119-
2120- real(dp)
2121- . A, BETA, D(2), DADD, DECUDD, DENMIN,
2122- . DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
2123- . DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
2124- . DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
2125- . DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
2126- . EC, ECUNIF, EX, EXUNIF,
2127- . F, F1, F2, F3, F4, FC, FX, FOUTHD,
2128- . GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
2129- . H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
2130- . T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
2131-
2132-C Lower bounds of density and its gradient to avoid divisions by zero
2133- PARAMETER ( DENMIN = 1.D-12 )
2134- PARAMETER ( GDMIN = 1.D-12 )
2135-
2136-C Fix some numerical parameters
2137- PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
2138- . THD=1.D0/3.D0, THRHLF=1.5D0,
2139- . TWO=2.D0, TWOTHD=2.D0/3.D0 )
2140-
2141-C Fix some more numerical constants
2142- PI = 4 * ATAN(1.D0)
2143- BETA = 0.066725D0
2144- GAMMA = (1 - LOG(TWO)) / PI**2
2145- MU = BETA * PI**2 / 3
2146- KAPPA = 0.804D0
2147-
2148-C Translate density and its gradient to new variables
2149- IF (nspin .EQ. 1) THEN
2150- D(1) = HALF * Dens(1)
2151- D(2) = D(1)
2152- DT = MAX( DENMIN, Dens(1) )
2153- DO 10 IX = 1,3
2154- GD(IX,1) = HALF * GDens(IX,1)
2155- GD(IX,2) = GD(IX,1)
2156- GDT(IX) = GDens(IX,1)
2157- 10 CONTINUE
2158- ELSE
2159- D(1) = Dens(1)
2160- D(2) = Dens(2)
2161- DT = MAX( DENMIN, Dens(1)+Dens(2) )
2162- DO 20 IX = 1,3
2163- GD(IX,1) = GDens(IX,1)
2164- GD(IX,2) = GDens(IX,2)
2165- GDT(IX) = GDens(IX,1) + GDens(IX,2)
2166- 20 CONTINUE
2167- ENDIF
2168- GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
2169- GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
2170- GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
2171- GDMT = MAX( GDMIN, GDMT )
2172-
2173-C Find local correlation energy and potential
2174- CALL PW92C( 2, D, ECUNIF, VCUNIF )
2175-
2176-C Find total correlation energy
2177- RS = ( 3 / (4*PI*DT) )**THD
2178- KF = (3 * PI**2 * DT)**THD
2179- KS = SQRT( 4 * KF / PI )
2180- ZETA = ( D(1) - D(2) ) / DT
2181- ZETA = MAX( -1.D0+DENMIN, ZETA )
2182- ZETA = MIN( 1.D0-DENMIN, ZETA )
2183- PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
2184- T = GDMT / (2 * PHI * KS * DT)
2185- F1 = ECUNIF / GAMMA / PHI**3
2186- F2 = EXP(-F1)
2187- A = BETA / GAMMA / (F2-1)
2188- F3 = T**2 + A * T**4
2189- F4 = BETA/GAMMA * F3 / (1 + A*F3)
2190- H = GAMMA * PHI**3 * LOG( 1 + F4 )
2191- FC = ECUNIF + H
2192-
2193-C Find correlation energy derivatives
2194- DRSDD = - (THD * RS / DT)
2195- DKFDD = THD * KF / DT
2196- DKSDD = HALF * KS * DKFDD / KF
2197- DZDD(1) = 1 / DT - ZETA / DT
2198- DZDD(2) = - (1 / DT) - ZETA / DT
2199- DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
2200- DO 40 IS = 1,2
2201- DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
2202- DPDD = DPDZ * DZDD(IS)
2203- DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
2204- DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
2205- DF2DD = (- F2) * DF1DD
2206- DADD = (- A) * DF2DD / (F2-1)
2207- DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
2208- DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
2209- DHDD = 3 * H * DPDD / PHI
2210- DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
2211- DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
2212-
2213- DO 30 IX = 1,3
2214- DTDGD = (T / GDMT) * GDT(IX) / GDMT
2215- DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
2216- DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
2217- DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
2218- DFCDGD(IX,IS) = DT * DHDGD
2219- 30 CONTINUE
2220- 40 CONTINUE
2221-
2222-C Find exchange energy and potential
2223- FX = 0
2224- DO 60 IS = 1,2
2225- DS(IS) = MAX( DENMIN, 2 * D(IS) )
2226- GDMS = MAX( GDMIN, 2 * GDM(IS) )
2227- KFS = (3 * PI**2 * DS(IS))**THD
2228- S = GDMS / (2 * KFS * DS(IS))
2229- F1 = 1 + MU * S**2 / KAPPA
2230- F = 1 + KAPPA - KAPPA / F1
2231-c
2232-c Note nspin=1 in call to exchng...
2233-c
2234- CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
2235- FX = FX + DS(IS) * EXUNIF * F
2236-
2237- DKFDD = THD * KFS / DS(IS)
2238- DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
2239- DF1DD = 2 * (F1-1) * DSDD / S
2240- DFDD = KAPPA * DF1DD / F1**2
2241- DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
2242-
2243- DO 50 IX = 1,3
2244- GDS = 2 * GD(IX,IS)
2245- DSDGD = (S / GDMS) * GDS / GDMS
2246- DF1DGD = 2 * MU * S * DSDGD / KAPPA
2247- DFDGD = KAPPA * DF1DGD / F1**2
2248- DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
2249- 50 CONTINUE
2250- 60 CONTINUE
2251- FX = HALF * FX / DT
2252-
2253-C Set output arguments
2254- EX = FX
2255- EC = FC
2256- DO 90 IS = 1,nspin
2257- DEXDD(IS) = DFXDD(IS)
2258- DECDD(IS) = DFCDD(IS)
2259- DO 80 IX = 1,3
2260- DEXDGD(IX,IS) = DFXDGD(IX,IS)
2261- DECDGD(IX,IS) = DFCDGD(IX,IS)
2262- 80 CONTINUE
2263- 90 CONTINUE
2264-
2265- END
2266-
2267- SUBROUTINE REVPBEXC( IREL, nspin, Dens, GDens,
2268- . EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
2269-
2270-C *********************************************************************
2271-C Implements revPBE: revised Perdew-Burke-Ernzerhof GGA.
2272-C Ref: Y. Zhang & W. Yang, Phys. Rev. Lett. 80, 890 (1998).
2273-C Written by E. Artacho in January 2006 by modifying the PBE routine of
2274-C L.C.Balbas and J.M.Soler. December 1996. Version 0.5.
2275-C ******** INPUT ******************************************************
2276-C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
2277-C INTEGER nspin : Number of spin polarizations (1 or 2)
2278-C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
2279-C spin electron density (if nspin=2)
2280-C REAL*8 GDens(3,nspin) : Total or spin density gradient
2281-C ******** OUTPUT *****************************************************
2282-C REAL*8 EX : Exchange energy density
2283-C REAL*8 EC : Correlation energy density
2284-C REAL*8 DEXDD(nspin) : Partial derivative
2285-C d(DensTot*Ex)/dDens(ispin),
2286-C where DensTot = Sum_ispin( Dens(ispin) )
2287-C For a constant density, this is the
2288-C exchange potential
2289-C REAL*8 DECDD(nspin) : Partial derivative
2290-C d(DensTot*Ec)/dDens(ispin),
2291-C where DensTot = Sum_ispin( Dens(ispin) )
2292-C For a constant density, this is the
2293-C correlation potential
2294-C REAL*8 DEXDGD(3,nspin): Partial derivative
2295-C d(DensTot*Ex)/d(GradDens(i,ispin))
2296-C REAL*8 DECDGD(3,nspin): Partial derivative
2297-C d(DensTot*Ec)/d(GradDens(i,ispin))
2298-C ********* UNITS ****************************************************
2299-C Lengths in Bohr
2300-C Densities in electrons per Bohr**3
2301-C Energies in Hartrees
2302-C Gradient vectors in cartesian coordinates
2303-C ********* ROUTINES CALLED ******************************************
2304-C EXCHNG, PW92C
2305-C ********************************************************************
2306-
2307- use precision, only : dp
2308-
2309- implicit none
2310- INTEGER IREL, nspin
2311- real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
2312- . DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
2313-
2314-C Internal variables
2315- INTEGER
2316- . IS, IX
2317-
2318- real(dp)
2319- . A, BETA, D(2), DADD, DECUDD, DENMIN,
2320- . DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
2321- . DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
2322- . DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
2323- . DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
2324- . EC, ECUNIF, EX, EXUNIF,
2325- . F, F1, F2, F3, F4, FC, FX, FOUTHD,
2326- . GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
2327- . H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
2328- . T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
2329-
2330-C Lower bounds of density and its gradient to avoid divisions by zero
2331- PARAMETER ( DENMIN = 1.D-12 )
2332- PARAMETER ( GDMIN = 1.D-12 )
2333-
2334-C Fix some numerical parameters
2335- PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
2336- . THD=1.D0/3.D0, THRHLF=1.5D0,
2337- . TWO=2.D0, TWOTHD=2.D0/3.D0 )
2338-
2339-C Fix some more numerical constants
2340- PI = 4 * ATAN(1.D0)
2341- BETA = 0.066725D0
2342- GAMMA = (1 - LOG(TWO)) / PI**2
2343- MU = BETA * PI**2 / 3
2344-cea The only modification w.r.t. PBE in this following line.
2345- KAPPA = 1.245D0
2346-
2347-C Translate density and its gradient to new variables
2348- IF (nspin .EQ. 1) THEN
2349- D(1) = HALF * Dens(1)
2350- D(2) = D(1)
2351- DT = MAX( DENMIN, Dens(1) )
2352- DO 10 IX = 1,3
2353- GD(IX,1) = HALF * GDens(IX,1)
2354- GD(IX,2) = GD(IX,1)
2355- GDT(IX) = GDens(IX,1)
2356- 10 CONTINUE
2357- ELSE
2358- D(1) = Dens(1)
2359- D(2) = Dens(2)
2360- DT = MAX( DENMIN, Dens(1)+Dens(2) )
2361- DO 20 IX = 1,3
2362- GD(IX,1) = GDens(IX,1)
2363- GD(IX,2) = GDens(IX,2)
2364- GDT(IX) = GDens(IX,1) + GDens(IX,2)
2365- 20 CONTINUE
2366- ENDIF
2367- GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
2368- GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
2369- GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
2370- GDMT = MAX( GDMIN, GDMT )
2371-
2372-C Find local correlation energy and potential
2373- CALL PW92C( 2, D, ECUNIF, VCUNIF )
2374-
2375-C Find total correlation energy
2376- RS = ( 3 / (4*PI*DT) )**THD
2377- KF = (3 * PI**2 * DT)**THD
2378- KS = SQRT( 4 * KF / PI )
2379- ZETA = ( D(1) - D(2) ) / DT
2380- ZETA = MAX( -1.D0+DENMIN, ZETA )
2381- ZETA = MIN( 1.D0-DENMIN, ZETA )
2382- PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
2383- T = GDMT / (2 * PHI * KS * DT)
2384- F1 = ECUNIF / GAMMA / PHI**3
2385- F2 = EXP(-F1)
2386- A = BETA / GAMMA / (F2-1)
2387- F3 = T**2 + A * T**4
2388- F4 = BETA/GAMMA * F3 / (1 + A*F3)
2389- H = GAMMA * PHI**3 * LOG( 1 + F4 )
2390- FC = ECUNIF + H
2391-
2392-C Find correlation energy derivatives
2393- DRSDD = - (THD * RS / DT)
2394- DKFDD = THD * KF / DT
2395- DKSDD = HALF * KS * DKFDD / KF
2396- DZDD(1) = 1 / DT - ZETA / DT
2397- DZDD(2) = - (1 / DT) - ZETA / DT
2398- DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
2399- DO 40 IS = 1,2
2400- DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
2401- DPDD = DPDZ * DZDD(IS)
2402- DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
2403- DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
2404- DF2DD = (- F2) * DF1DD
2405- DADD = (- A) * DF2DD / (F2-1)
2406- DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
2407- DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
2408- DHDD = 3 * H * DPDD / PHI
2409- DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
2410- DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
2411-
2412- DO 30 IX = 1,3
2413- DTDGD = (T / GDMT) * GDT(IX) / GDMT
2414- DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
2415- DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
2416- DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
2417- DFCDGD(IX,IS) = DT * DHDGD
2418- 30 CONTINUE
2419- 40 CONTINUE
2420-
2421-C Find exchange energy and potential
2422- FX = 0
2423- DO 60 IS = 1,2
2424- DS(IS) = MAX( DENMIN, 2 * D(IS) )
2425- GDMS = MAX( GDMIN, 2 * GDM(IS) )
2426- KFS = (3 * PI**2 * DS(IS))**THD
2427- S = GDMS / (2 * KFS * DS(IS))
2428- F1 = 1 + MU * S**2 / KAPPA
2429- F = 1 + KAPPA - KAPPA / F1
2430-c
2431-c Note nspin=1 in call to exchng...
2432-c
2433- CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
2434- FX = FX + DS(IS) * EXUNIF * F
2435-
2436- DKFDD = THD * KFS / DS(IS)
2437- DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
2438- DF1DD = 2 * (F1-1) * DSDD / S
2439- DFDD = KAPPA * DF1DD / F1**2
2440- DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
2441-
2442- DO 50 IX = 1,3
2443- GDS = 2 * GD(IX,IS)
2444- DSDGD = (S / GDMS) * GDS / GDMS
2445- DF1DGD = 2 * MU * S * DSDGD / KAPPA
2446- DFDGD = KAPPA * DF1DGD / F1**2
2447- DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
2448- 50 CONTINUE
2449- 60 CONTINUE
2450- FX = HALF * FX / DT
2451-
2452-C Set output arguments
2453- EX = FX
2454- EC = FC
2455- DO 90 IS = 1,nspin
2456- DEXDD(IS) = DFXDD(IS)
2457- DECDD(IS) = DFCDD(IS)
2458- DO 80 IX = 1,3
2459- DEXDGD(IX,IS) = DFXDGD(IX,IS)
2460- DECDGD(IX,IS) = DFCDGD(IX,IS)
2461- 80 CONTINUE
2462- 90 CONTINUE
2463-
2464- END
2465-
2466-
2467- SUBROUTINE PW91XC( IREL, nspin, Dens, GDens,
2468- . EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
2469-
2470-C *********************************************************************
2471-C Implements Perdew-Wang91 Generalized-Gradient-Approximation.
2472-C Ref: JCP 100, 1290 (1994)
2473-C Written by J.L. Martins August 2000
2474-C ******** INPUT ******************************************************
2475-C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
2476-C INTEGER nspin : Number of spin polarizations (1 or 2)
2477-C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
2478-C spin electron density (if nspin=2)
2479-C REAL*8 GDens(3,nspin) : Total or spin density gradient
2480-C ******** OUTPUT *****************************************************
2481-C REAL*8 EX : Exchange energy density
2482-C REAL*8 EC : Correlation energy density
2483-C REAL*8 DEXDD(nspin) : Partial derivative
2484-C d(DensTot*Ex)/dDens(ispin),
2485-C where DensTot = Sum_ispin( Dens(ispin) )
2486-C For a constant density, this is the
2487-C exchange potential
2488-C REAL*8 DECDD(nspin) : Partial derivative
2489-C d(DensTot*Ec)/dDens(ispin),
2490-C where DensTot = Sum_ispin( Dens(ispin) )
2491-C For a constant density, this is the
2492-C correlation potential
2493-C REAL*8 DEXDGD(3,nspin): Partial derivative
2494-C d(DensTot*Ex)/d(GradDens(i,ispin))
2495-C REAL*8 DECDGD(3,nspin): Partial derivative
2496-C d(DensTot*Ec)/d(GradDens(i,ispin))
2497-C ********* UNITS ****************************************************
2498-C Lengths in Bohr
2499-C Densities in electrons per Bohr**3
2500-C Energies in Hartrees
2501-C Gradient vectors in cartesian coordinates
2502-C ********* ROUTINES CALLED ******************************************
2503-C EXCHNG, PW92C
2504-C ********************************************************************
2505-
2506- use precision, only : dp
2507-
2508- implicit none
2509- INTEGER IREL, nspin
2510- real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
2511- . DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
2512-
2513-C Internal variables
2514- INTEGER
2515- . IS, IX
2516- real(dp)
2517- . A, BETA, D(2), DADD, DECUDD, DENMIN,
2518- . DF1DD, DF2DD, DF3DD, DF4DD, DF3DGD, DF4DGD,
2519- . DFCDD(2), DFCDGD(3,2), DFDGD, DFXDD(2), DFXDGD(3,2),
2520- . DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
2521- . DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
2522- . EC, ECUNIF, EX, EXUNIF,
2523- . F, F1, F2, F3, F4, FC, FX, FOUTHD,
2524- . GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
2525- . H, HALF, KF, KFS, KS, PHI, PI, RS, S,
2526- . T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
2527-
2528- real(dp) F5, F6, F7, F8, ASINHS
2529- real(dp) DF5DD,DF6DD,DF7DD,DF8DD
2530- real(dp) DF1DS, DF2DS, DF3DS, DFDS, DF7DGD
2531-
2532-C Lower bounds of density and its gradient to avoid divisions by zero
2533- PARAMETER ( DENMIN = 1.D-12 )
2534- PARAMETER ( GDMIN = 1.D-12 )
2535-
2536-C Fix some numerical parameters
2537- PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
2538- . THD=1.D0/3.D0, THRHLF=1.5D0,
2539- . TWO=2.D0, TWOTHD=2.D0/3.D0 )
2540-
2541-C Fix some more numerical constants
2542- PI = 4.0_dp * ATAN(1.0_dp)
2543- BETA = 15.75592_dp * 0.004235_dp
2544- GAMMA = BETA**2 / (2.0_dp * 0.09_dp)
2545-
2546-C Translate density and its gradient to new variables
2547- IF (nspin .EQ. 1) THEN
2548- D(1) = HALF * Dens(1)
2549- D(2) = D(1)
2550- DT = MAX( DENMIN, Dens(1) )
2551- DO 10 IX = 1,3
2552- GD(IX,1) = HALF * GDens(IX,1)
2553- GD(IX,2) = GD(IX,1)
2554- GDT(IX) = GDens(IX,1)
2555- 10 CONTINUE
2556- ELSE
2557- D(1) = Dens(1)
2558- D(2) = Dens(2)
2559- DT = MAX( DENMIN, Dens(1)+Dens(2) )
2560- DO 20 IX = 1,3
2561- GD(IX,1) = GDens(IX,1)
2562- GD(IX,2) = GDens(IX,2)
2563- GDT(IX) = GDens(IX,1) + GDens(IX,2)
2564- 20 CONTINUE
2565- ENDIF
2566- GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
2567- GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
2568- GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
2569- GDMT = MAX( GDMIN, GDMT )
2570-
2571-C Find local correlation energy and potential
2572- CALL PW92C( 2, D, ECUNIF, VCUNIF )
2573-
2574-C Find total correlation energy
2575- RS = ( 3 / (4*PI*DT) )**THD
2576- KF = (3 * PI**2 * DT)**THD
2577- KS = SQRT( 4 * KF / PI )
2578- S = GDMT / (2 * KF * DT)
2579- T = GDMT / (2 * KS * DT)
2580- ZETA = ( D(1) - D(2) ) / DT
2581- ZETA = MAX( -1.D0+DENMIN, ZETA )
2582- ZETA = MIN( 1.D0-DENMIN, ZETA )
2583- PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
2584- F1 = ECUNIF / GAMMA / PHI**3
2585- F2 = EXP(-F1)
2586- A = BETA / GAMMA / (F2-1)
2587- F3 = T**2 + A * T**4
2588- F4 = BETA/GAMMA * F3 / (1 + A*F3)
2589- F5 = 0.002568D0 + 0.023266D0*RS + 7.389D-6*RS**2
2590- F6 = 1.0D0 + 8.723D0*RS + 0.472D0*RS**2 + 0.07389D0*RS**3
2591- F7 = EXP(-100.0D0 * S**2 * PHI**4)
2592- F8 = 15.75592D0*(0.001667212D0 + F5/F6 -0.004235D0 +
2593- . 3.0D0*0.001667212D0/7.0D0)
2594- H = GAMMA * PHI**3 * LOG( 1 + F4 ) + F8 * T**2 * F7
2595- FC = ECUNIF + H
2596-
2597-C Find correlation energy derivatives
2598- DRSDD = - THD * RS / DT
2599- DKFDD = THD * KF / DT
2600- DKSDD = HALF * KS * DKFDD / KF
2601- DZDD(1) = 1 / DT - ZETA / DT
2602- DZDD(2) = - 1 / DT - ZETA / DT
2603- DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
2604- DO 40 IS = 1,2
2605- DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
2606- DPDD = DPDZ * DZDD(IS)
2607- DTDD = - T * ( DPDD/PHI + DKSDD/KS + 1/DT )
2608- DSDD = - S * ( DPDD/PHI + DKFDD/KF + 1/DT )
2609- DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
2610- DF2DD = - F2 * DF1DD
2611- DADD = - A * DF2DD / (F2-1)
2612- DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
2613- DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
2614- DF5DD = (0.023266D0 + 2.0D0*7.389D-6*RS)*DRSDD
2615- DF6DD = (8.723D0 + 2.0D0*0.472D0*RS
2616- . + 3.0D0*0.07389D0*RS**2)*DRSDD
2617- DF7DD = -200.0D0 * S * PHI**4 * DSDD * F7
2618- . -100.0D0 * S**2 * 4.0D0* PHI**3 * DPDD * F7
2619- DF8DD = 15.75592D0 * DF5DD/F6 - 15.75592D0*F5*DF6DD / F6**2
2620- DHDD = 3 * H * DPDD / PHI
2621- DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
2622- DHDD = DHDD + DF8DD * T**2 * F7
2623- DHDD = DHDD + F8 * 2*T*DTDD *F7
2624- DHDD = DHDD + F8 * T**2 * DF7DD
2625-
2626- DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
2627- DO 30 IX = 1,3
2628- DTDGD = (T / GDMT) * GDT(IX) / GDMT
2629- DSDGD = (S / GDMT) * GDT(IX) / GDMT
2630- DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
2631- DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
2632- DF7DGD = -200.0D0 * S * PHI**4 * DSDGD * F7
2633- DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
2634- DHDGD = DHDGD + F8 * 2*T*DTDGD *F7 + F8 * T**2 *DF7DGD
2635- DFCDGD(IX,IS) = DT * DHDGD
2636- 30 CONTINUE
2637- 40 CONTINUE
2638-
2639-C Find exchange energy and potential
2640- FX = 0
2641- DO 60 IS = 1,2
2642- DS(1) = MAX( DENMIN, 2 * D(IS) )
2643- GDMS = MAX( GDMIN, 2 * GDM(IS) )
2644- KFS = (3 * PI**2 * DS(1))**THD
2645- S = GDMS / (2 * KFS * DS(1))
2646- F4 = SQRT(1.0D0 + (7.7956D0*S)**2)
2647- ASINHS = LOG(7.7956D0*S + F4)
2648- F1 = 1.0D0 + 0.19645D0 * S * ASINHS
2649- F2 = 0.2743D0 - 0.15084D0*EXP(-100.0D0*S*S)
2650- F3 = 1.0D0 / (F1 + 0.004D0 * S*S*S*S)
2651- F = (F1 + F2 * S*S ) * F3
2652- .
2653- CALL EXCHNG( IREL, 1, DS, EXUNIF, VXUNIF )
2654- FX = FX + DS(1) * EXUNIF * F
2655-
2656- DKFDD = THD * KFS / DS(1)
2657- DSDD = S * ( -DKFDD/KFS - 1/DS(1) )
2658- DF1DS = 0.19645D0 * ASINHS +
2659- . 0.19645D0 * S * 7.7956D0 / F4
2660- DF2DS = 0.15084D0*200.0D0*S*EXP(-100.0D0*S*S)
2661- DF3DS = - F3*F3 * (DF1DS + 4.0D0*0.004D0 * S*S*S)
2662- DFDS = DF1DS * F3 + DF2DS * S*S * F3 + 2.0D0 * S * F2 * F3
2663- . + (F1 + F2 * S*S ) * DF3DS
2664- DFXDD(IS) = VXUNIF(1) * F + DS(1) * EXUNIF * DFDS * DSDD
2665-
2666- DO 50 IX = 1,3
2667- GDS = 2 * GD(IX,IS)
2668- DSDGD = (S / GDMS) * GDS / GDMS
2669- DFDGD = DFDS * DSDGD
2670- DFXDGD(IX,IS) = DS(1) * EXUNIF * DFDGD
2671- 50 CONTINUE
2672- 60 CONTINUE
2673- FX = HALF * FX / DT
2674-
2675-C Set output arguments
2676- EX = FX
2677- EC = FC
2678- DO 90 IS = 1,nspin
2679- DEXDD(IS) = DFXDD(IS)
2680- DECDD(IS) = DFCDD(IS)
2681- DO 80 IX = 1,3
2682- DEXDGD(IX,IS) = DFXDGD(IX,IS)
2683- DECDGD(IX,IS) = DFCDGD(IX,IS)
2684- 80 CONTINUE
2685- 90 CONTINUE
2686-
2687- END
2688-
2689-
2690-
2691- SUBROUTINE PW92C( nspin, Dens, EC, VC )
2692-
2693-C ********************************************************************
2694-C Implements the Perdew-Wang'92 local correlation (beyond RPA).
2695-C Ref: J.P.Perdew & Y.Wang, PRB, 45, 13244 (1992)
2696-C Written by L.C.Balbas and J.M.Soler. Dec'96. Version 0.5.
2697-C ********* INPUT ****************************************************
2698-C INTEGER nspin : Number of spin polarizations (1 or 2)
2699-C REAL*8 Dens(nspin) : Local (spin) density
2700-C ********* OUTPUT ***************************************************
2701-C REAL*8 EC : Correlation energy density
2702-C REAL*8 VC(nspin) : Correlation (spin) potential
2703-C ********* UNITS ****************************************************
2704-C Densities in electrons per Bohr**3
2705-C Energies in Hartrees
2706-C ********* ROUTINES CALLED ******************************************
2707-C None
2708-C ********************************************************************
2709-
2710- use precision, only : dp
2711-
2712-C Next line is nonstandard but may be supressed
2713- implicit none
2714-
2715-C Argument types and dimensions
2716- INTEGER nspin
2717- real(dp) Dens(nspin), EC, VC(nspin)
2718-
2719-C Internal variable declarations
2720- INTEGER IG
2721- real(dp) A(0:2), ALPHA1(0:2), B, BETA(0:2,4), C,
2722- . DBDRS, DECDD(2), DECDRS, DECDZ, DENMIN, DFDZ,
2723- . DGDRS(0:2), DCDRS, DRSDD, DTOT, DZDD(2),
2724- . F, FPP0, FOUTHD, G(0:2), HALF, ONE,
2725- . P(0:2), PI, RS, THD, THRHLF, ZETA
2726-
2727-C Add tiny numbers to avoid numerical errors
2728- PARAMETER ( DENMIN = 1.D-12 )
2729- PARAMETER ( ONE = 1.D0 + 1.D-12 )
2730-
2731-C Fix some numerical constants
2732- PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
2733- . THD=1.D0/3.D0, THRHLF=1.5D0 )
2734-
2735-C Parameters from Table I of Perdew & Wang, PRB, 45, 13244 (92)
2736- DATA P / 1.00d0, 1.00d0, 1.00d0 /
2737- DATA A / 0.031091d0, 0.015545d0, 0.016887d0 /
2738- DATA ALPHA1 / 0.21370d0, 0.20548d0, 0.11125d0 /
2739- DATA BETA / 7.5957d0, 14.1189d0, 10.357d0,
2740- . 3.5876d0, 6.1977d0, 3.6231d0,
2741- . 1.6382d0, 3.3662d0, 0.88026d0,
2742- . 0.49294d0, 0.62517d0, 0.49671d0 /
2743-
2744-C Find rs and zeta
2745- PI = 4 * ATAN(1.D0)
2746- IF (nspin .EQ. 1) THEN
2747- DTOT = MAX( DENMIN, Dens(1) )
2748- ZETA = 0
2749- RS = ( 3 / (4*PI*DTOT) )**THD
2750-C Find derivatives dRs/dDens and dZeta/dDens
2751- DRSDD = (- RS) / DTOT / 3
2752- DZDD(1) = 0
2753- ELSE
2754- DTOT = MAX( DENMIN, Dens(1)+Dens(2) )
2755- ZETA = ( Dens(1) - Dens(2) ) / DTOT
2756- RS = ( 3 / (4*PI*DTOT) )**THD
2757- DRSDD = (- RS) / DTOT / 3
2758- DZDD(1) = (ONE - ZETA) / DTOT
2759- DZDD(2) = - (ONE + ZETA) / DTOT
2760- ENDIF
2761-
2762-C Find eps_c(rs,0)=G(0), eps_c(rs,1)=G(1) and -alpha_c(rs)=G(2)
2763-C using eq.(10) of cited reference (Perdew & Wang, PRB, 45, 13244 (92))
2764- DO 20 IG = 0,2
2765- B = BETA(IG,1) * RS**HALF +
2766- . BETA(IG,2) * RS +
2767- . BETA(IG,3) * RS**THRHLF +
2768- . BETA(IG,4) * RS**(P(IG)+1)
2769- DBDRS = BETA(IG,1) * HALF / RS**HALF +
2770- . BETA(IG,2) +
2771- . BETA(IG,3) * THRHLF * RS**HALF +
2772- . BETA(IG,4) * (P(IG)+1) * RS**P(IG)
2773- C = 1 + 1 / (2 * A(IG) * B)
2774- DCDRS = - ( (C-1) * DBDRS / B )
2775- G(IG) = (- 2) * A(IG) * ( 1 + ALPHA1(IG)*RS ) * LOG(C)
2776- DGDRS(IG) = (- 2) *A(IG) * ( ALPHA1(IG) * LOG(C) +
2777- . (1+ALPHA1(IG)*RS) * DCDRS / C )
2778- 20 CONTINUE
2779-
2780-C Find f''(0) and f(zeta) from eq.(9)
2781- C = 1 / (2**FOUTHD - 2)
2782- FPP0 = 8 * C / 9
2783- F = ( (ONE+ZETA)**FOUTHD + (ONE-ZETA)**FOUTHD - 2 ) * C
2784- DFDZ = FOUTHD * ( (ONE+ZETA)**THD - (ONE-ZETA)**THD ) * C
2785-
2786-C Find eps_c(rs,zeta) from eq.(8)
2787- EC = G(0) - G(2) * F / FPP0 * (ONE-ZETA**4) +
2788- . (G(1)-G(0)) * F * ZETA**4
2789- DECDRS = DGDRS(0) - DGDRS(2) * F / FPP0 * (ONE-ZETA**4) +
2790- . (DGDRS(1)-DGDRS(0)) * F * ZETA**4
2791- DECDZ = (- G(2)) / FPP0 * ( DFDZ*(ONE-ZETA**4) - F*4*ZETA**3 ) +
2792- . (G(1)-G(0)) * ( DFDZ*ZETA**4 + F*4*ZETA**3 )
2793-
2794-C Find correlation potential
2795- IF (nspin .EQ. 1) THEN
2796- DECDD(1) = DECDRS * DRSDD
2797- VC(1) = EC + DTOT * DECDD(1)
2798- ELSE
2799- DECDD(1) = DECDRS * DRSDD + DECDZ * DZDD(1)
2800- DECDD(2) = DECDRS * DRSDD + DECDZ * DZDD(2)
2801- VC(1) = EC + DTOT * DECDD(1)
2802- VC(2) = EC + DTOT * DECDD(2)
2803- ENDIF
2804-
2805- END
2806-
2807-
2808-
2809- SUBROUTINE PW92XC( IREL, nspin, Dens, EPSX, EPSC, VX, VC )
2810-
2811-C ********************************************************************
2812-C Implements the Perdew-Wang'92 LDA/LSD exchange correlation
2813-C Ref: J.P.Perdew & Y.Wang, PRB, 45, 13244 (1992)
2814-C Written by L.C.Balbas and J.M.Soler. Dec'96. Version 0.5.
2815-C ********* INPUT ****************************************************
2816-C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
2817-C INTEGER nspin : Number of spin polarizations (1 or 2)
2818-C REAL*8 Dens(nspin) : Local (spin) density
2819-C ********* OUTPUT ***************************************************
2820-C REAL*8 EPSX : Exchange energy density
2821-C REAL*8 EPSC : Correlation energy density
2822-C REAL*8 VX(nspin) : Exchange (spin) potential
2823-C REAL*8 VC(nspin) : Correlation (spin) potential
2824-C ********* UNITS ****************************************************
2825-C Densities in electrons per Bohr**3
2826-C Energies in Hartrees
2827-C ********* ROUTINES CALLED ******************************************
2828-C EXCHNG, PW92C
2829-C ********************************************************************
2830-
2831- use precision, only : dp
2832-
2833- implicit none
2834- INTEGER IREL, nspin
2835- real(dp) Dens(nspin), EPSX, EPSC, VC(nspin), VX(nspin)
2836-
2837- CALL EXCHNG( IREL, nspin, Dens, EPSX, VX )
2838- CALL PW92C( nspin, Dens, EPSC, VC )
2839- END
2840-
2841-
2842-
2843- SUBROUTINE PZXC( IREL, NSP, DS, EX, EC, VX, VC, DVXDN, DVCDN )
2844-
2845-C *****************************************************************
2846-C Perdew-Zunger parameterization of Ceperley-Alder exchange and
2847-C correlation. Ref: Perdew & Zunger, Phys. Rev. B 23 5075 (1981).
2848-C Adapted by J.M.Soler from routine velect of Froyen's
2849-C pseudopotential generation program. Madrid, Jan'97.
2850-C **** Input *****************************************************
2851-C INTEGER IREL : relativistic-exchange switch (0=no, 1=yes)
2852-C INTEGER NSP : spin-polarizations (1=>unpolarized, 2=>polarized)
2853-C REAL*8 DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
2854-C **** Output *****************************************************
2855-C REAL*8 EX : exchange energy density
2856-C REAL*8 EC : correlation energy density
2857-C REAL*8 VX(NSP) : (spin-dependent) exchange potential
2858-C REAL*8 VC(NSP) : (spin-dependent) correlation potential
2859-C REAL*8 DVXDN(NSP,NSP): Derivative of the exchange potential
2860-C respect the charge density,
2861-C Dvx(spin1)/Dn(spin2)
2862-C REAL*8 DVCDN(NSP,NSP): Derivative of the correlation potential
2863-C respect the charge density,
2864-C Dvc(spin1)/Dn(spin2)
2865-C **** Units *******************************************************
2866-C Densities in electrons/Bohr**3
2867-C Energies in Hartrees
2868-C *****************************************************************
2869-
2870- use precision, only: dp
2871-
2872- implicit none
2873-
2874- integer :: nsp, irel, isp1, isp2, isp
2875- real(dp) :: DS(NSP), VX(NSP), VC(NSP),
2876- . DVXDN(NSP,NSP), DVCDN(NSP,NSP)
2877- real(dp), parameter ::
2878- $ ZERO=0.D0,ONE=1.D0,PFIVE=.5D0,OPF=1.5D0,PNN=.99D0,
2879- $ PTHREE=0.3D0,PSEVF=0.75D0,C0504=0.0504D0,
2880- $ C0254=0.0254D0,C014=0.014D0,C0406=0.0406D0,
2881- $ C15P9=15.9D0,C0666=0.0666D0,C11P4=11.4D0,
2882- $ C045=0.045D0,C7P8=7.8D0,C88=0.88D0,C20P59=20.592D0,
2883- $ C3P52=3.52D0,C0311=0.0311D0,C0014=0.0014D0,
2884- $ C0538=0.0538D0,C0096=0.0096D0,C096=0.096D0,
2885- $ C0622=0.0622D0,C004=0.004D0,C0232=0.0232D0,
2886- $ C1686=0.1686D0,C1P398=1.3981D0,C2611=0.2611D0,
2887- $ C2846=0.2846D0,C1P053=1.0529D0,C3334=0.3334D0
2888-
2889-C Ceperly-Alder 'ca' constants. Internal energies in Rydbergs.
2890- real(dp), parameter ::
2891- $ CON1=1.D0/6, CON2=0.008D0/3, CON3=0.3502D0/3,
2892- $ CON4=0.0504D0/3, CON5=0.0028D0/3, CON6=0.1925D0/3,
2893- $ CON7=0.0206D0/3, CON8=9.7867D0/6, CON9=1.0444D0/3,
2894- $ CON10=7.3703D0/6, CON11=1.3336D0/3
2895-
2896-C X-alpha parameter:
2897- real(dp), PARAMETER :: ALP = 2.D0 / 3.D0
2898-
2899-C Other variables converted into parameters by J.M.Soler
2900- real(dp), parameter ::
2901- $ TINY = 1.D-6 ,
2902- $ PI = 3.14159265358979312_dp,
2903- $ TWO = 2.0D0,
2904- $ HALF = 0.5D0,
2905- $ TRD = 1.D0 / 3.D0,
2906- $ FTRD = 4.D0 / 3.D0,
2907- $ TFTM = 0.51984209978974638D0,
2908- $ A0 = 0.52106176119784808D0,
2909- $ CRS = 0.620350490899400087D0,
2910- $ CXP = (- 3.D0) * ALP / (PI*A0),
2911- $ CXF = 1.25992104989487319D0
2912-
2913- real(dp) :: d1, d2, d, z, fz, fzp
2914- real(dp) :: ex, ec, dfzpdn, rs, vxp, exp_var
2915- real(dp) :: beta, sb, alb, vxf, exf, dvxpdn
2916- real(dp) :: dvxfdn, sqrs, te, be, ecp, vcp
2917- real(dp) :: dtedn, be2, dbedn, dvcpdn, decpdn
2918- real(dp) :: ecf, vcf, dvcfdn, decfdn, rslog
2919-
2920-
2921-C Find density and polarization
2922- IF (NSP .EQ. 2) THEN
2923- D1 = MAX(DS(1),ZERO)
2924- D2 = MAX(DS(2),ZERO)
2925- D = D1 + D2
2926- IF (D .LE. ZERO) THEN
2927- EX = ZERO
2928- EC = ZERO
2929- VX(1) = ZERO
2930- VX(2) = ZERO
2931- VC(1) = ZERO
2932- VC(2) = ZERO
2933- RETURN
2934- ENDIF
2935-c
2936-c Robustness enhancement by Jose Soler (August 2002)
2937-c
2938- Z = (D1 - D2) / D
2939- IF (Z .LE. -ONE) THEN
2940- FZ = (TWO**FTRD-TWO)/TFTM
2941- FZP = -FTRD*TWO**TRD/TFTM
2942- DFZPDN = FTRD*TRD*TWO**(-ALP)/TFTM
2943- ELSEIF (Z .GE. ONE) THEN
2944- FZ = (TWO**FTRD-TWO)/TFTM
2945- FZP = FTRD*TWO**TRD/TFTM
2946- DFZPDN = FTRD*TRD*TWO**(-ALP)/TFTM
2947- ELSE
2948- FZ = ((ONE+Z)**FTRD+(ONE-Z)**FTRD-TWO)/TFTM
2949- FZP = FTRD*((ONE+Z)**TRD-(ONE-Z)**TRD)/TFTM
2950- DFZPDN = FTRD*TRD*((ONE+Z)**(-ALP) + (ONE-Z)**(-ALP))/TFTM
2951- ENDIF
2952- ELSE
2953- D = DS(1)
2954- IF (D .LE. ZERO) THEN
2955- EX = ZERO
2956- EC = ZERO
2957- VX(1) = ZERO
2958- VC(1) = ZERO
2959- RETURN
2960- ENDIF
2961- Z = ZERO
2962- FZ = ZERO
2963- FZP = ZERO
2964- ENDIF
2965- RS = CRS / D**TRD
2966-
2967-C Exchange
2968- VXP = CXP / RS
2969- EXP_VAR = 0.75D0 * VXP
2970- IF (IREL .EQ. 1) THEN
2971- BETA = C014/RS
2972- IF (BETA .LT. TINY) THEN
2973- SB = ONE + HALF*BETA**2
2974- ALB = BETA
2975- ELSE
2976- SB = SQRT(1+BETA*BETA)
2977- ALB = LOG(BETA+SB)
2978- ENDIF
2979- VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
2980- EXP_VAR = EXP_VAR *(ONE-OPF*((BETA*SB-ALB)/BETA**2)**2)
2981- ENDIF
2982- VXF = CXF * VXP
2983- EXF = CXF * EXP_VAR
2984- DVXPDN = TRD * VXP / D
2985- DVXFDN = TRD * VXF / D
2986-
2987-C Correlation
2988- IF (RS .GT. ONE) THEN
2989- SQRS=SQRT(RS)
2990- TE = ONE+CON10*SQRS+CON11*RS
2991- BE = ONE+C1P053*SQRS+C3334*RS
2992- ECP = -(C2846/BE)
2993- VCP = ECP*TE/BE
2994- DTEDN = ((CON10 * SQRS *HALF) + CON11 * RS)*(-TRD/D)
2995- BE2 = BE * BE
2996- DBEDN = ((C1P053 * SQRS *HALF) + C3334 * RS)*(-TRD/D)
2997- DVCPDN = -(C2846/BE2)*(DTEDN - 2.0D0 * TE * DBEDN/BE)
2998- DECPDN = (C2846/BE2)*DBEDN
2999- TE = ONE+CON8*SQRS+CON9*RS
3000- BE = ONE+C1P398*SQRS+C2611*RS
3001- ECF = -(C1686/BE)
3002- VCF = ECF*TE/BE
3003- DTEDN = ((CON8 * SQRS * HALF) + CON9 * RS)*(-TRD/D)
3004- BE2 = BE * BE
3005- DBEDN = ((C1P398 * SQRS * HALF) + C2611 * RS)*(-TRD/D)
3006- DVCFDN = -(C1686/BE2)*(DTEDN - 2.0D0 * TE * DBEDN/BE)
3007- DECFDN = (C1686/BE2)*DBEDN
3008- ELSE
3009- RSLOG=LOG(RS)
3010- ECP=(C0622+C004*RS)*RSLOG-C096-C0232*RS
3011- VCP=(C0622+CON2*RS)*RSLOG-CON3-CON4*RS
3012- DVCPDN = (CON2*RS*RSLOG + (CON2-CON4)*RS + C0622)*(-TRD/D)
3013- DECPDN = (C004*RS*RSLOG + (C004-C0232)*RS + C0622)*(-TRD/D)
3014- ECF=(C0311+C0014*RS)*RSLOG-C0538-C0096*RS
3015- VCF=(C0311+CON5*RS)*RSLOG-CON6-CON7*RS
3016- DVCFDN = (CON5*RS*RSLOG + (CON5-CON7)*RS + C0311)*(-TRD/D)
3017- DECFDN = (C0014*RS*RSLOG + (C0014-C0096)*RS + C0311)*(-TRD/D)
3018- ENDIF
3019-
3020- ISP1 = 1
3021- ISP2 = 2
3022-
3023-C Find up and down potentials
3024- IF (NSP .EQ. 2) THEN
3025- EX = EXP_VAR + FZ*(EXF-EXP_VAR)
3026- EC = ECP + FZ*(ECF-ECP)
3027- VX(1) = VXP + FZ*(VXF-VXP) + (ONE-Z)*FZP*(EXF-EXP_VAR)
3028- VX(2) = VXP + FZ*(VXF-VXP) - (ONE+Z)*FZP*(EXF-EXP_VAR)
3029- VC(1) = VCP + FZ*(VCF-VCP) + (ONE-Z)*FZP*(ECF-ECP)
3030- VC(2) = VCP + FZ*(VCF-VCP) - (ONE+Z)*FZP*(ECF-ECP)
3031-
3032-C Derivatives of exchange potential respect the density
3033-
3034- DVXDN(ISP1,ISP1) =
3035- . DVXPDN
3036- . + FZP*(VXF-VXP-EXF+EXP_VAR)*( 2.D0*D2/(D*D) )
3037- . + FZ*(DVXFDN-DVXPDN)+(1-Z)*FZP*(VXF-VXP)/(4.D0*D)
3038- . + (1-Z)*DFZPDN*(EXF-EXP_VAR)*( 2.D0*D2/(D*D) )
3039- DVXDN(ISP1,ISP2) =
3040- . DVXPDN
3041- . + FZP*(VXF-VXP-EXF+EXP_VAR)*(-2.D0*D1/(D*D) )
3042- . + FZ*(DVXFDN-DVXPDN)+(1-Z)*FZP*(VXF-VXP)/(4.D0*D)
3043- . + (1-Z)*DFZPDN*(EXF-EXP_VAR)*( -2.D0*D1/(D*D) )
3044- DVXDN(ISP2,ISP1) =
3045- . DVXPDN
3046- . + FZP*(VXF-VXP-EXF+EXP_VAR)*( 2.D0*D2/(D*D) )
3047- . + FZ*(DVXFDN-DVXPDN)-(1+Z)*FZP*(VXF-VXP)/(4.D0*D)
3048- . - (1+Z)*DFZPDN*(EXF-EXP_VAR)*( 2.D0*D2/(D*D) )
3049- DVXDN(ISP2,ISP2) =
3050- . DVXPDN
3051- . + FZP*(VXF-VXP-EXF+EXP_VAR)*(-2.D0*D1/(D*D) )
3052- . + FZ*(DVXFDN-DVXPDN)-(1+Z)*FZP*(VXF-VXP)/(4.D0*D)
3053- . - (1+Z)*DFZPDN*(EXF-EXP_VAR)*( -2.D0*D1/(D*D) )
3054-
3055-C Derivatives of correlation potential respect the density
3056-
3057- DVCDN(ISP1,ISP1) =
3058- . DVCPDN
3059- . + FZP*(VCF-VCP-ECF+ECP)*( 2.D0*D2/(D*D) )
3060- . + FZ*(DVCFDN-DVCPDN)+ (1-Z)*FZP*(DECFDN-DECPDN)
3061- . + (1-Z)*DFZPDN*(ECF-ECP)*( 2.D0*D2/(D*D) )
3062- DVCDN(ISP1,ISP2) =
3063- . DVCPDN
3064- . + FZP*(VCF-VCP-ECF+ECP)*(-2.D0*D1/(D*D) )
3065- . + FZ*(DVCFDN-DVCPDN)+ (1-Z)*FZP*(DECFDN-DECPDN)
3066- . + (1-Z)*DFZPDN*(ECF-ECP)*( -2.D0*D1/(D*D) )
3067- DVCDN(ISP2,ISP1) =
3068- . DVCPDN
3069- . + FZP*(VCF-VCP-ECF+ECP)*( 2.D0*D2/(D*D) )
3070- . + FZ*(DVCFDN-DVCPDN)- (1+Z)*FZP*(DECFDN-DECPDN)
3071- . - (1+Z)*DFZPDN*(ECF-ECP)*( 2.D0*D2/(D*D) )
3072- DVCDN(ISP2,ISP2) =
3073- . DVCPDN
3074- . + FZP*(VCF-VCP-ECF+ECP)*(-2.D0*D1/(D*D) )
3075- . + FZ*(DVCFDN-DVCPDN)- (1+Z)*FZP*(DECFDN-DECPDN)
3076- . - (1+Z)*DFZPDN*(ECF-ECP)*( -2.D0*D1/(D*D) )
3077-
3078- ELSE
3079- EX = EXP_VAR
3080- EC = ECP
3081- VX(1) = VXP
3082- VC(1) = VCP
3083- DVXDN(1,1) = DVXPDN
3084- DVCDN(1,1) = DVCPDN
3085- ENDIF
3086-
3087-C Change from Rydbergs to Hartrees
3088- EX = HALF * EX
3089- EC = HALF * EC
3090- DO 10 ISP = 1,NSP
3091- VX(ISP) = HALF * VX(ISP)
3092- VC(ISP) = HALF * VC(ISP)
3093- DO 5 ISP2 = 1,NSP
3094- DVXDN(ISP,ISP2) = HALF * DVXDN(ISP,ISP2)
3095- DVCDN(ISP,ISP2) = HALF * DVCDN(ISP,ISP2)
3096- 5 CONTINUE
3097- 10 CONTINUE
3098- END
3099-
3100- subroutine blypxc(nspin,dens,gdens,EX,EC,
3101- . dEXdd,dECdd,dEXdgd,dECdgd)
3102-c ***************************************************************
3103-c Implements Becke gradient exchange functional (A.D.
3104-c Becke, Phys. Rev. A 38, 3098 (1988)) and Lee, Yang, Parr
3105-c correlation functional (C. Lee, W. Yang, R.G. Parr, Phys. Rev. B
3106-c 37, 785 (1988)), as modificated by Miehlich,Savin,Stoll and Preuss,
3107-c Chem. Phys. Lett. 157,200 (1989). See also Johnson, Gill and Pople,
3108-c J. Chem. Phys. 98, 5612 (1993). Some errors were detected in this
3109-c last paper, so not all of the expressions correspond exactly to those
3110-c implemented here.
3111-c Written by Maider Machado. July 1998.
3112-c **************** INPUT ********************************************
3113-c integer nspin : Number of spin polarizations (1 or 2)
3114-c real*8 dens(nspin) : Total electron density (if nspin=1) or
3115-c spin electron density (if nspin=2)
3116-c real*8 gdens(3,nspin) : Total or spin density gradient
3117-c ******** OUTPUT *****************************************************
3118-c real*8 ex : Exchange energy density
3119-c real*8 ec : Correlation energy density
3120-c real*8 dexdd(nspin) : Partial derivative
3121-c d(DensTot*Ex)/dDens(ispin),
3122-c where DensTot = Sum_ispin( Dens(ispin) )
3123-c For a constant density, this is the
3124-c exchange potential
3125-c real*8 decdd(nspin) : Partial derivative
3126-c d(DensTot*Ec)/dDens(ispin),
3127-c where DensTot = Sum_ispin( Dens(ispin) )
3128-c For a constant density, this is the
3129-c correlation potential
3130-c real*8 dexdgd(3,nspin): Partial derivative
3131-c d(DensTot*Ex)/d(GradDens(i,ispin))
3132-c real*8 decdgd(3,nspin): Partial derivative
3133-c d(DensTot*Ec)/d(GradDens(i,ispin))
3134-c ********* UNITS ****************************************************
3135-c Lengths in Bohr
3136-c Densities in electrons per Bohr**3
3137-c Energies in Hartrees
3138-c Gradient vectors in cartesian coordinates
3139-c ********************************************************************
3140-
3141- use precision, only : dp
3142-
3143- implicit none
3144-
3145- integer nspin
3146- real(dp) dens(nspin), gdens(3,nspin), EX, EC,
3147- . dEXdd(nspin), dECdd(nspin), dEXdgd(3,nspin),
3148- . dECdgd(3,nspin)
3149-
3150-c Internal variables
3151- integer is,ix
3152- real(dp) pi, beta, thd, tthd, thrhlf, half, fothd,
3153- . d(2),gd(3,2),dmin, ash,gdm(2),denmin,dt,
3154- . g(2),x(2),a,b,c,dd,onzthd,gdmin,
3155- . ga, gb, gc,becke,dbecgd(3,2),
3156- . dgdx(2), dgdxa, dgdxb, dgdxc,dgdxd,dbecdd(2),
3157- . den,omega, domega, delta, ddelta,cf,
3158- . gam11, gam12, gam22, LYPa, LYPb1,
3159- . LYPb2,dLYP11,dLYP12,dLYP22,LYP,
3160- . dd1g11,dd1g12,dd1g22,dd2g12,dd2g11,dd2g22,
3161- . dLYPdd(2),dg11dd(3,2),dg22dd(3,2),
3162- . dLYPgd(3,2)
3163-
3164-c Lower bounds of density and its gradient to avoid divisions by zero
3165- parameter ( denmin=1.d-8 )
3166- parameter (gdmin=1.d-8)
3167- parameter (dmin=1.d-5)
3168-
3169-c Fix some numerical parameters
3170- parameter ( thd = 1.d0/3.d0, tthd=2.d0/3.d0 )
3171- parameter ( thrhlf=1.5d0, half=0.5d0,
3172- . fothd=4.d0/3.d0, onzthd=11.d0/3.d0)
3173-
3174-c Empirical parameter for Becke exchange functional (a.u.)
3175- parameter(beta= 0.0042d0)
3176-
3177-c Constants for LYP functional (a.u.)
3178- parameter(a=0.04918d0, b=0.132d0, c=0.2533d0, dd=0.349d0)
3179-
3180- pi= 4*atan(1.d0)
3181-
3182-
3183-c Translate density and its gradient to new variables
3184- if (nspin .eq. 1) then
3185- d(1) = half * dens(1)
3186- d(1) = max(denmin,d(1))
3187- d(2) = d(1)
3188- dt = max( denmin, dens(1) )
3189- do ix = 1,3
3190- gd(ix,1) = half * gdens(ix,1)
3191- gd(ix,2) = gd(ix,1)
3192- enddo
3193- else
3194- d(1) = dens(1)
3195- d(2) = dens(2)
3196- do is=1,2
3197- d(is) = max (denmin,d(is))
3198- enddo
3199- dt = max( denmin, dens(1)+dens(2) )
3200- do ix = 1,3
3201- gd(ix,1) = gdens(ix,1)
3202- gd(ix,2) = gdens(ix,2)
3203- enddo
3204- endif
3205-
3206- gdm(1) = sqrt( gd(1,1)**2 + gd(2,1)**2 + gd(3,1)**2 )
3207- gdm(2) = sqrt( gd(1,2)**2 + gd(2,2)**2 + gd(3,2)**2 )
3208-
3209- do is=1,2
3210- gdm(is)= max(gdm(is),gdmin)
3211- enddo
3212-
3213-c Find Becke exchange energy
3214- ga = -thrhlf*(3.d0/4.d0/pi)**thd
3215- do is=1,2
3216- if(d(is).lt.dmin) then
3217- g(is)=ga
3218- else
3219- x(is) = gdm(is)/d(is)**fothd
3220- gb = beta*x(is)**2
3221- ash=log(x(is)+sqrt(x(is)**2+1))
3222- gc = 1+6*beta*x(is)*ash
3223- g(is) = ga-gb/gc
3224- endif
3225- enddo
3226-
3227-c Density of energy
3228- becke=(g(1)*d(1)**fothd+g(2)*d(2)**fothd)/dt
3229-
3230-
3231-c Exchange energy derivatives
3232- do is=1,2
3233- if(d(is).lt.dmin)then
3234- dbecdd(is)=0.
3235- do ix=1,3
3236- dbecgd(ix,is)=0.
3237- enddo
3238- else
3239- dgdxa=6*beta**2*x(is)**2
3240- ash=log(x(is)+sqrt(x(is)**2+1))
3241- dgdxb=x(is)/sqrt(x(is)**2+1)-ash
3242- dgdxc=-2*beta*x(is)
3243- dgdxd=(1+6*beta*x(is)*ash)**2
3244- dgdx(is)=(dgdxa*dgdxb+dgdxc)/dgdxd
3245- dbecdd(is)=fothd*d(is)**thd*(g(is)-x(is)*dgdx(is))
3246- do ix=1,3
3247- dbecgd(ix,is)=d(is)**(-fothd)*dgdx(is)*gd(ix,is)/x(is)
3248- enddo
3249- endif
3250- enddo
3251-
3252-c Lee-Yang-Parr correlation energy
3253- den=1+dd*dt**(-thd)
3254- omega=dt**(-onzthd)*exp(-c*dt**(-thd))/den
3255- delta=c*dt**(-thd)+dd*dt**(-thd)/den
3256- cf=3.*(3*pi**2)**tthd/10.
3257- gam11=gdm(1)**2
3258- gam12=gd(1,1)*gd(1,2)+gd(2,1)*gd(2,2)+gd(3,1)*gd(3,2)
3259- gam22=gdm(2)**2
3260- LYPa=-4*a*d(1)*d(2)/(den*dt)
3261- LYPb1=2**onzthd*cf*a*b*omega*d(1)*d(2)
3262- LYPb2=d(1)**(8./3.)+d(2)**(8./3.)
3263- dLYP11=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)
3264- .*d(1)/dt)-d(2)**2)
3265- dLYP12=-a*b*omega*(d(1)*d(2)/9.*(47.-7.*delta)
3266- .-fothd*dt**2)
3267- dLYP22=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)*
3268- .d(2)/dt)-d(1)**2)
3269-
3270-c Density of energy
3271- LYP=(LYPa-LYPb1*LYPb2+dLYP11*gam11+dLYP12*gam12
3272- .+dLYP22*gam22)/dt
3273-
3274-c Correlation energy derivatives
3275- domega=-thd*dt**(-fothd)*omega*(11.*dt**thd-c-dd/den)
3276- ddelta=thd*(dd**2*dt**(-5./3.)/den**2-delta/dt)
3277-
3278-c Second derivatives with respect to the density
3279- dd1g11=domega/omega*dLYP11-a*b*omega*(d(2)/9.*
3280- . (1.-3.*delta-2*(delta-11.)*d(1)/dt)-d(1)*d(2)/9.*
3281- . ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2))
3282-
3283- dd1g12=domega/omega*dLYP12-a*b*omega*(d(2)/9.*
3284- . (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
3285-
3286- dd1g22=domega/omega*dLYP22-a*b*omega*(1./9.*d(2)
3287- . *(1.-3.*delta-(delta-11.)*d(2)/dt)-d(1)*d(2)/9.*
3288- . ((3.+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2)-2*d(1))
3289-
3290-
3291- dd2g22=domega/omega*dLYP22-a*b*omega*(d(1)/9.*
3292- . (1.-3.*delta-2*(delta-11.)*d(2)/dt)-d(1)*d(2)/9.*
3293- . ((3+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2))
3294-
3295-
3296- dd2g12=domega/omega*dLYP12-a*b*omega*(d(1)/9.*
3297- . (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt)
3298-
3299- dd2g11=domega/omega*dLYP11-a*b*omega*(1./9.*d(1)
3300- . *(1.-3.*delta-(delta-11.)*d(1)/dt)-d(1)*d(2)/9.*
3301- . ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2)-2*d(2))
3302-
3303-
3304- dLYPdd(1)=-4*a/den*d(1)*d(2)/dt*
3305- . (thd*dd*dt**(-fothd)/den
3306- . +1./d(1)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)*
3307- . (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(2)*(onzthd*
3308- . d(1)**(8./3.)+d(2)**(8./3.)))+dd1g11*gam11+
3309- . dd1g12*gam12+dd1g22*gam22
3310-
3311-
3312- dLYPdd(2)=-4*a/den*d(1)*d(2)/dt*(thd*dd*dt**(-fothd)/den
3313- . +1./d(2)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)*
3314- . (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(1)*(onzthd*
3315- . d(2)**(8./3.)+d(1)**(8./3.)))+dd2g22*gam22+
3316- . dd2g12*gam12+dd2g11*gam11
3317-
3318-
3319-c second derivatives with respect to the density gradient
3320-
3321- do is=1,2
3322- do ix=1,3
3323- dg11dd(ix,is)=2*gd(ix,is)
3324- dg22dd(ix,is)=2*gd(ix,is)
3325- enddo
3326- enddo
3327- do ix=1,3
3328- dLYPgd(ix,1)=dLYP11*dg11dd(ix,1)+dLYP12*gd(ix,2)
3329- dLYPgd(ix,2)=dLYP22*dg22dd(ix,2)+dLYP12*gd(ix,1)
3330- enddo
3331-
3332-
3333- EX=becke
3334- EC=LYP
3335- do is=1,nspin
3336- dEXdd(is)=dbecdd(is)
3337- dECdd(is)=dLYPdd(is)
3338- do ix=1,3
3339- dEXdgd(ix,is)=dbecgd(ix,is)
3340- dECdgd(ix,is)=dLYPgd(ix,is)
3341- enddo
3342- enddo
3343- end
3344-
3345- SUBROUTINE RPBEXC( IREL, nspin, Dens, GDens,
3346- . EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
3347-
3348-C *********************************************************************
3349-C Implements Hammer's RPBE Generalized-Gradient-Approximation (GGA).
3350-C A revision of PBE (Perdew-Burke-Ernzerhof)
3351-C Ref: Hammer, Hansen & Norskov, PRB 59, 7413 (1999) and
3352-C J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996)
3353-C
3354-C Written by M.V. Fernandez-Serra. March 2004. On the PBE routine of
3355-C L.C.Balbas and J.M.Soler. December 1996. Version 0.5.
3356-C ******** INPUT ******************************************************
3357-C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
3358-C INTEGER nspin : Number of spin polarizations (1 or 2)
3359-C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
3360-C spin electron density (if nspin=2)
3361-C REAL*8 GDens(3,nspin) : Total or spin density gradient
3362-C ******** OUTPUT *****************************************************
3363-C REAL*8 EX : Exchange energy density
3364-C REAL*8 EC : Correlation energy density
3365-C REAL*8 DEXDD(nspin) : Partial derivative
3366-C d(DensTot*Ex)/dDens(ispin),
3367-C where DensTot = Sum_ispin( Dens(ispin) )
3368-C For a constant density, this is the
3369-C exchange potential
3370-C REAL*8 DECDD(nspin) : Partial derivative
3371-C d(DensTot*Ec)/dDens(ispin),
3372-C where DensTot = Sum_ispin( Dens(ispin) )
3373-C For a constant density, this is the
3374-C correlation potential
3375-C REAL*8 DEXDGD(3,nspin): Partial derivative
3376-C d(DensTot*Ex)/d(GradDens(i,ispin))
3377-C REAL*8 DECDGD(3,nspin): Partial derivative
3378-C d(DensTot*Ec)/d(GradDens(i,ispin))
3379-C ********* UNITS ****************************************************
3380-C Lengths in Bohr
3381-C Densities in electrons per Bohr**3
3382-C Energies in Hartrees
3383-C Gradient vectors in cartesian coordinates
3384-C ********* ROUTINES CALLED ******************************************
3385-C EXCHNG, PW92C
3386-C ********************************************************************
3387-
3388- use precision, only : dp
3389-
3390- implicit none
3391- INTEGER IREL, nspin
3392- real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
3393- . DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
3394-
3395-C Internal variables
3396- INTEGER
3397- . IS, IX
3398-
3399- real(dp)
3400- . A, BETA, D(2), DADD, DECUDD, DENMIN,
3401- . DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
3402- . DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
3403- . DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
3404- . DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
3405- . EC, ECUNIF, EX, EXUNIF,
3406- . F, F1, F2, F3, F4, FC, FX, FOUTHD,
3407- . GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
3408- . H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
3409- . T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
3410-
3411-C Lower bounds of density and its gradient to avoid divisions by zero
3412- PARAMETER ( DENMIN = 1.D-12 )
3413- PARAMETER ( GDMIN = 1.D-12 )
3414-
3415-C Fix some numerical parameters
3416- PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
3417- . THD=1.D0/3.D0, THRHLF=1.5D0,
3418- . TWO=2.D0, TWOTHD=2.D0/3.D0 )
3419-
3420-C Fix some more numerical constants
3421- PI = 4 * ATAN(1.D0)
3422- BETA = 0.066725D0
3423- GAMMA = (1 - LOG(TWO)) / PI**2
3424- MU = BETA * PI**2 / 3
3425- KAPPA = 0.804D0
3426-
3427-C Translate density and its gradient to new variables
3428- IF (nspin .EQ. 1) THEN
3429- D(1) = HALF * Dens(1)
3430- D(2) = D(1)
3431- DT = MAX( DENMIN, Dens(1) )
3432- DO 10 IX = 1,3
3433- GD(IX,1) = HALF * GDens(IX,1)
3434- GD(IX,2) = GD(IX,1)
3435- GDT(IX) = GDens(IX,1)
3436- 10 CONTINUE
3437- ELSE
3438- D(1) = Dens(1)
3439- D(2) = Dens(2)
3440- DT = MAX( DENMIN, Dens(1)+Dens(2) )
3441- DO 20 IX = 1,3
3442- GD(IX,1) = GDens(IX,1)
3443- GD(IX,2) = GDens(IX,2)
3444- GDT(IX) = GDens(IX,1) + GDens(IX,2)
3445- 20 CONTINUE
3446- ENDIF
3447- GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
3448- GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
3449- GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
3450- GDMT = MAX( GDMIN, GDMT )
3451-
3452-C Find local correlation energy and potential
3453- CALL PW92C( 2, D, ECUNIF, VCUNIF )
3454-
3455-C Find total correlation energy
3456- RS = ( 3 / (4*PI*DT) )**THD
3457- KF = (3 * PI**2 * DT)**THD
3458- KS = SQRT( 4 * KF / PI )
3459- ZETA = ( D(1) - D(2) ) / DT
3460- ZETA = MAX( -1.D0+DENMIN, ZETA )
3461- ZETA = MIN( 1.D0-DENMIN, ZETA )
3462- PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
3463- T = GDMT / (2 * PHI * KS * DT)
3464- F1 = ECUNIF / GAMMA / PHI**3
3465- F2 = EXP(-F1)
3466- A = BETA / GAMMA / (F2-1)
3467- F3 = T**2 + A * T**4
3468- F4 = BETA/GAMMA * F3 / (1 + A*F3)
3469- H = GAMMA * PHI**3 * LOG( 1 + F4 )
3470- FC = ECUNIF + H
3471-
3472-C Find correlation energy derivatives
3473- DRSDD = - (THD * RS / DT)
3474- DKFDD = THD * KF / DT
3475- DKSDD = HALF * KS * DKFDD / KF
3476- DZDD(1) = 1 / DT - ZETA / DT
3477- DZDD(2) = - (1 / DT) - ZETA / DT
3478- DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
3479- DO 40 IS = 1,2
3480- DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
3481- DPDD = DPDZ * DZDD(IS)
3482- DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
3483- DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
3484- DF2DD = (- F2) * DF1DD
3485- DADD = (- A) * DF2DD / (F2-1)
3486- DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
3487- DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
3488- DHDD = 3 * H * DPDD / PHI
3489- DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
3490- DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
3491-
3492- DO 30 IX = 1,3
3493- DTDGD = (T / GDMT) * GDT(IX) / GDMT
3494- DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
3495- DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
3496- DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
3497- DFCDGD(IX,IS) = DT * DHDGD
3498- 30 CONTINUE
3499- 40 CONTINUE
3500-
3501-C Find exchange energy and potential
3502- FX = 0
3503- DO 60 IS = 1,2
3504- DS(IS) = MAX( DENMIN, 2 * D(IS) )
3505- GDMS = MAX( GDMIN, 2 * GDM(IS) )
3506- KFS = (3 * PI**2 * DS(IS))**THD
3507- S = GDMS / (2 * KFS * DS(IS))
3508-cea Hammer's RPBE (Hammer, Hansen & Norskov PRB 59 7413 (99)
3509-cea F1 = DEXP( - MU * S**2 / KAPPA)
3510-cea F = 1 + KAPPA * (1 - F1)
3511-cea Following is standard PBE
3512-cea F1 = 1 + MU * S**2 / KAPPA
3513-cea F = 1 + KAPPA - KAPPA / F1
3514-cea (If revPBE Zhang & Yang, PRL 80,890(1998),change PBE's KAPPA to 1.245)
3515- F1 = DEXP( - MU * S**2 / KAPPA)
3516- F = 1 + KAPPA * (1 - F1)
3517-
3518-c Note nspin=1 in call to exchng...
3519-
3520- CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
3521- FX = FX + DS(IS) * EXUNIF * F
3522-
3523-cMVFS The derivatives of F also need to be changed for Hammer's RPBE.
3524-cMVFS DF1DD = 2 * F1 * DSDD * ( - MU * S / KAPPA)
3525-cMVFS DF1DGD= 2 * F1 * DSDGD * ( - MU * S / KAPPA)
3526-cMVFS DFDD = -1 * KAPPA * DF1DD
3527-cMVFS DFDGD = -1 * KAPPA * DFDGD
3528-
3529- DKFDD = THD * KFS / DS(IS)
3530- DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
3531-c DF1DD = 2 * (F1-1) * DSDD / S
3532-c DFDD = KAPPA * DF1DD / F1**2
3533- DF1DD = 2* F1 * DSDD * ( - MU * S / KAPPA)
3534- DFDD = -1 * KAPPA * DF1DD
3535- DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
3536-
3537- DO 50 IX = 1,3
3538- GDS = 2 * GD(IX,IS)
3539- DSDGD = (S / GDMS) * GDS / GDMS
3540-c DF1DGD = 2 * MU * S * DSDGD / KAPPA
3541-c DFDGD = KAPPA * DF1DGD / F1**2
3542- DF1DGD =2*F1 * DSDGD * ( - MU * S / KAPPA)
3543- DFDGD = -1 * KAPPA * DF1DGD
3544- DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
3545- 50 CONTINUE
3546- 60 CONTINUE
3547- FX = HALF * FX / DT
3548-
3549-C Set output arguments
3550- EX = FX
3551- EC = FC
3552- DO 90 IS = 1,nspin
3553- DEXDD(IS) = DFXDD(IS)
3554- DECDD(IS) = DFCDD(IS)
3555- DO 80 IX = 1,3
3556- DEXDGD(IX,IS) = DFXDGD(IX,IS)
3557- DECDGD(IX,IS) = DFCDGD(IX,IS)
3558- 80 CONTINUE
3559- 90 CONTINUE
3560-
3561- END
3562- SUBROUTINE WCXC( IREL, nspin, Dens, GDens,
3563- . EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
3564-
3565-C *********************************************************************
3566-C Implements Wu-Cohen Generalized-Gradient-Approximation.
3567-C Ref: Z. Wu and R. E. Cohen PRB 73, 235116 (2006)
3568-C Written by Marivi Fernandez-Serra, with contributions by
3569-C Julian Gale and Alberto Garcia,
3570-C over the PBEXC subroutine of L.C.Balbas and J.M.Soler.
3571-C September, 2006.
3572-C ******** INPUT ******************************************************
3573-C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
3574-C INTEGER nspin : Number of spin polarizations (1 or 2)
3575-C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
3576-C spin electron density (if nspin=2)
3577-C REAL*8 GDens(3,nspin) : Total or spin density gradient
3578-C ******** OUTPUT *****************************************************
3579-C REAL*8 EX : Exchange energy density
3580-C REAL*8 EC : Correlation energy density
3581-C REAL*8 DEXDD(nspin) : Partial derivative
3582-C d(DensTot*Ex)/dDens(ispin),
3583-C where DensTot = Sum_ispin( Dens(ispin) )
3584-C For a constant density, this is the
3585-C exchange potential
3586-C REAL*8 DECDD(nspin) : Partial derivative
3587-C d(DensTot*Ec)/dDens(ispin),
3588-C where DensTot = Sum_ispin( Dens(ispin) )
3589-C For a constant density, this is the
3590-C correlation potential
3591-C REAL*8 DEXDGD(3,nspin): Partial derivative
3592-C d(DensTot*Ex)/d(GradDens(i,ispin))
3593-C REAL*8 DECDGD(3,nspin): Partial derivative
3594-C d(DensTot*Ec)/d(GradDens(i,ispin))
3595-C ********* UNITS ****************************************************
3596-C Lengths in Bohr
3597-C Densities in electrons per Bohr**3
3598-C Energies in Hartrees
3599-C Gradient vectors in cartesian coordinates
3600-C ********* ROUTINES CALLED ******************************************
3601-C EXCHNG, PW92C
3602-C ********************************************************************
3603-
3604- use precision, only : dp
3605-
3606- implicit none
3607- INTEGER IREL, nspin
3608- real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
3609- . DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
3610-
3611-C Internal variables
3612- INTEGER
3613- . IS, IX
3614-
3615- real(dp)
3616- . A, BETA, D(2), DADD, DECUDD, DENMIN,
3617- . DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
3618- . DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
3619- . DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
3620- . DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
3621- . XWC, DXWCDS, CWC,
3622- . EC, ECUNIF, EX, EXUNIF,
3623- . F, F1, F2, F3, F4, FC, FX, FOUTHD,
3624- . GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
3625- . H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
3626- . TEN81,
3627- . T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
3628-
3629-C Lower bounds of density and its gradient to avoid divisions by zero
3630- PARAMETER ( DENMIN = 1.D-12 )
3631- PARAMETER ( GDMIN = 1.D-12 )
3632-
3633-C Fix some numerical parameters
3634- PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
3635- . THD=1.D0/3.D0, THRHLF=1.5D0,
3636- . TWO=2.D0, TWOTHD=2.D0/3.D0 )
3637- PARAMETER ( TEN81 = 10.0d0/81.0d0 )
3638-
3639-C Fix some more numerical constants
3640- PI = 4 * ATAN(1.D0)
3641- BETA = 0.066725D0
3642- GAMMA = (1 - LOG(TWO)) / PI**2
3643- MU = BETA * PI**2 / 3
3644- KAPPA = 0.804D0
3645- CWC = 0.0079325D0
3646-
3647-C Translate density and its gradient to new variables
3648- IF (nspin .EQ. 1) THEN
3649- D(1) = HALF * Dens(1)
3650- D(2) = D(1)
3651- DT = MAX( DENMIN, Dens(1) )
3652- DO 10 IX = 1,3
3653- GD(IX,1) = HALF * GDens(IX,1)
3654- GD(IX,2) = GD(IX,1)
3655- GDT(IX) = GDens(IX,1)
3656- 10 CONTINUE
3657- ELSE
3658- D(1) = Dens(1)
3659- D(2) = Dens(2)
3660- DT = MAX( DENMIN, Dens(1)+Dens(2) )
3661- DO 20 IX = 1,3
3662- GD(IX,1) = GDens(IX,1)
3663- GD(IX,2) = GDens(IX,2)
3664- GDT(IX) = GDens(IX,1) + GDens(IX,2)
3665- 20 CONTINUE
3666- ENDIF
3667- GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
3668- GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
3669- GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
3670- GDMT = MAX( GDMIN, GDMT )
3671-
3672-C Find local correlation energy and potential
3673- CALL PW92C( 2, D, ECUNIF, VCUNIF )
3674-
3675-C Find total correlation energy
3676- RS = ( 3 / (4*PI*DT) )**THD
3677- KF = (3 * PI**2 * DT)**THD
3678- KS = SQRT( 4 * KF / PI )
3679- ZETA = ( D(1) - D(2) ) / DT
3680- ZETA = MAX( -1.D0+DENMIN, ZETA )
3681- ZETA = MIN( 1.D0-DENMIN, ZETA )
3682- PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
3683- T = GDMT / (2 * PHI * KS * DT)
3684- F1 = ECUNIF / GAMMA / PHI**3
3685- F2 = EXP(-F1)
3686- A = BETA / GAMMA / (F2-1)
3687- F3 = T**2 + A * T**4
3688- F4 = BETA/GAMMA * F3 / (1 + A*F3)
3689- H = GAMMA * PHI**3 * LOG( 1 + F4 )
3690- FC = ECUNIF + H
3691-
3692-C Find correlation energy derivatives
3693- DRSDD = - (THD * RS / DT)
3694- DKFDD = THD * KF / DT
3695- DKSDD = HALF * KS * DKFDD / KF
3696- DZDD(1) = 1 / DT - ZETA / DT
3697- DZDD(2) = - (1 / DT) - ZETA / DT
3698- DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
3699- DO 40 IS = 1,2
3700- DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
3701- DPDD = DPDZ * DZDD(IS)
3702- DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
3703- DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
3704- DF2DD = (- F2) * DF1DD
3705- DADD = (- A) * DF2DD / (F2-1)
3706- DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
3707- DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
3708- DHDD = 3 * H * DPDD / PHI
3709- DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
3710- DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
3711-
3712- DO 30 IX = 1,3
3713- DTDGD = (T / GDMT) * GDT(IX) / GDMT
3714- DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
3715- DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
3716- DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
3717- DFCDGD(IX,IS) = DT * DHDGD
3718- 30 CONTINUE
3719- 40 CONTINUE
3720-
3721-C Find exchange energy and potential
3722- FX = 0
3723- DO 60 IS = 1,2
3724- DS(IS) = MAX( DENMIN, 2 * D(IS) )
3725- GDMS = MAX( GDMIN, 2 * GDM(IS) )
3726- KFS = (3 * PI**2 * DS(IS))**THD
3727- S = GDMS / (2 * KFS * DS(IS))
3728-c
3729-c For PBE:
3730-c
3731-c x = MU * S**2
3732-c dxds = 2*MU*S
3733-c
3734-c Wu-Cohen form:
3735-c
3736- XWC= TEN81 * s**2 + (MU- TEN81) *
3737- . S**2 * exp(-S**2) + log(1+ CWC * S**4)
3738- DXWCDS = 2 * TEN81 * S + (MU - TEN81) * exp(-S**2) *
3739- . 2*S * (1 - S*S) + 4 * CWC * S**3 / (1 + CWC * S**4)
3740-c-------------------
3741-
3742- F1 = 1 + XWC / KAPPA
3743- F = 1 + KAPPA - KAPPA / F1
3744-c
3745-c Note nspin=1 in call to exchng...
3746-c
3747- CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
3748- FX = FX + DS(IS) * EXUNIF * F
3749-
3750- DKFDD = THD * KFS / DS(IS)
3751- DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
3752- DF1DD = DXWCDS * DSDD / KAPPA
3753- DFDD = KAPPA * DF1DD / F1**2
3754- DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
3755-
3756- DO 50 IX = 1,3
3757- GDS = 2 * GD(IX,IS)
3758- DSDGD = (S / GDMS) * GDS / GDMS
3759- DF1DGD = DXWCDS * DSDGD / KAPPA
3760- DFDGD = KAPPA * DF1DGD / F1**2
3761- DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
3762- 50 CONTINUE
3763- 60 CONTINUE
3764- FX = HALF * FX / DT
3765-
3766-C Set output arguments
3767- EX = FX
3768- EC = FC
3769- DO 90 IS = 1,nspin
3770- DEXDD(IS) = DFXDD(IS)
3771- DECDD(IS) = DFCDD(IS)
3772- DO 80 IX = 1,3
3773- DEXDGD(IX,IS) = DFXDGD(IX,IS)
3774- DECDGD(IX,IS) = DFCDGD(IX,IS)
3775- 80 CONTINUE
3776- 90 CONTINUE
3777-
3778- END
3779-
3780- SUBROUTINE PBESOLXC( IREL, nspin, Dens, GDens,
3781- . EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
3782-
3783-C *********************************************************************
3784-C Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation.
3785-C with the revised parameters for solids (PBEsol).
3786-C Ref: J.P.Perdew et al, PRL 100, 136406 (2008)
3787-C Written by L.C.Balbas and J.M.Soler for PBE. December 1996.
3788-C Modified by J.D. Gale for PBEsol. May 2009.
3789-C ******** INPUT ******************************************************
3790-C INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes)
3791-C INTEGER nspin : Number of spin polarizations (1 or 2)
3792-C REAL*8 Dens(nspin) : Total electron density (if nspin=1) or
3793-C spin electron density (if nspin=2)
3794-C REAL*8 GDens(3,nspin) : Total or spin density gradient
3795-C ******** OUTPUT *****************************************************
3796-C REAL*8 EX : Exchange energy density
3797-C REAL*8 EC : Correlation energy density
3798-C REAL*8 DEXDD(nspin) : Partial derivative
3799-C d(DensTot*Ex)/dDens(ispin),
3800-C where DensTot = Sum_ispin( Dens(ispin) )
3801-C For a constant density, this is the
3802-C exchange potential
3803-C REAL*8 DECDD(nspin) : Partial derivative
3804-C d(DensTot*Ec)/dDens(ispin),
3805-C where DensTot = Sum_ispin( Dens(ispin) )
3806-C For a constant density, this is the
3807-C correlation potential
3808-C REAL*8 DEXDGD(3,nspin): Partial derivative
3809-C d(DensTot*Ex)/d(GradDens(i,ispin))
3810-C REAL*8 DECDGD(3,nspin): Partial derivative
3811-C d(DensTot*Ec)/d(GradDens(i,ispin))
3812-C ********* UNITS ****************************************************
3813-C Lengths in Bohr
3814-C Densities in electrons per Bohr**3
3815-C Energies in Hartrees
3816-C Gradient vectors in cartesian coordinates
3817-C ********* ROUTINES CALLED ******************************************
3818-C EXCHNG, PW92C
3819-C ********************************************************************
3820-
3821- use precision, only : dp
3822-
3823- implicit none
3824- INTEGER IREL, nspin
3825- real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
3826- . DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin)
3827-
3828-C Internal variables
3829- INTEGER
3830- . IS, IX
3831-
3832- real(dp)
3833- . A, BETA, D(2), DADD, DECUDD, DENMIN,
3834- . DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
3835- . DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
3836- . DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
3837- . DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
3838- . EC, ECUNIF, EX, EXUNIF,
3839- . F, F1, F2, F3, F4, FC, FX, FOUTHD,
3840- . GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
3841- . H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
3842- . T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA
3843-
3844-C Lower bounds of density and its gradient to avoid divisions by zero
3845- PARAMETER ( DENMIN = 1.D-12 )
3846- PARAMETER ( GDMIN = 1.D-12 )
3847-
3848-C Fix some numerical parameters
3849- PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
3850- . THD=1.D0/3.D0, THRHLF=1.5D0,
3851- . TWO=2.D0, TWOTHD=2.D0/3.D0 )
3852-
3853-C Fix some more numerical constants
3854- PI = 4 * ATAN(1.D0)
3855- BETA = 0.046d0
3856- GAMMA = (1 - LOG(TWO)) / PI**2
3857- MU = 10.0d0/81.0d0
3858- KAPPA = 0.804D0
3859-
3860-C Translate density and its gradient to new variables
3861- IF (nspin .EQ. 1) THEN
3862- D(1) = HALF * Dens(1)
3863- D(2) = D(1)
3864- DT = MAX( DENMIN, Dens(1) )
3865- DO 10 IX = 1,3
3866- GD(IX,1) = HALF * GDens(IX,1)
3867- GD(IX,2) = GD(IX,1)
3868- GDT(IX) = GDens(IX,1)
3869- 10 CONTINUE
3870- ELSE
3871- D(1) = Dens(1)
3872- D(2) = Dens(2)
3873- DT = MAX( DENMIN, Dens(1)+Dens(2) )
3874- DO 20 IX = 1,3
3875- GD(IX,1) = GDens(IX,1)
3876- GD(IX,2) = GDens(IX,2)
3877- GDT(IX) = GDens(IX,1) + GDens(IX,2)
3878- 20 CONTINUE
3879- ENDIF
3880- GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
3881- GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
3882- GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 )
3883- GDMT = MAX( GDMIN, GDMT )
3884-
3885-C Find local correlation energy and potential
3886- CALL PW92C( 2, D, ECUNIF, VCUNIF )
3887-
3888-C Find total correlation energy
3889- RS = ( 3 / (4*PI*DT) )**THD
3890- KF = (3 * PI**2 * DT)**THD
3891- KS = SQRT( 4 * KF / PI )
3892- ZETA = ( D(1) - D(2) ) / DT
3893- ZETA = MAX( -1.D0+DENMIN, ZETA )
3894- ZETA = MIN( 1.D0-DENMIN, ZETA )
3895- PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
3896- T = GDMT / (2 * PHI * KS * DT)
3897- F1 = ECUNIF / GAMMA / PHI**3
3898- F2 = EXP(-F1)
3899- A = BETA / GAMMA / (F2-1)
3900- F3 = T**2 + A * T**4
3901- F4 = BETA/GAMMA * F3 / (1 + A*F3)
3902- H = GAMMA * PHI**3 * LOG( 1 + F4 )
3903- FC = ECUNIF + H
3904-
3905-C Find correlation energy derivatives
3906- DRSDD = - (THD * RS / DT)
3907- DKFDD = THD * KF / DT
3908- DKSDD = HALF * KS * DKFDD / KF
3909- DZDD(1) = 1 / DT - ZETA / DT
3910- DZDD(2) = - (1 / DT) - ZETA / DT
3911- DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
3912- DO 40 IS = 1,2
3913- DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
3914- DPDD = DPDZ * DZDD(IS)
3915- DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
3916- DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
3917- DF2DD = (- F2) * DF1DD
3918- DADD = (- A) * DF2DD / (F2-1)
3919- DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
3920- DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )
3921- DHDD = 3 * H * DPDD / PHI
3922- DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4)
3923- DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
3924-
3925- DO 30 IX = 1,3
3926- DTDGD = (T / GDMT) * GDT(IX) / GDMT
3927- DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 )
3928- DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )
3929- DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4)
3930- DFCDGD(IX,IS) = DT * DHDGD
3931- 30 CONTINUE
3932- 40 CONTINUE
3933-
3934-C Find exchange energy and potential
3935- FX = 0
3936- DO 60 IS = 1,2
3937- DS(IS) = MAX( DENMIN, 2 * D(IS) )
3938- GDMS = MAX( GDMIN, 2 * GDM(IS) )
3939- KFS = (3 * PI**2 * DS(IS))**THD
3940- S = GDMS / (2 * KFS * DS(IS))
3941- F1 = 1 + MU * S**2 / KAPPA
3942- F = 1 + KAPPA - KAPPA / F1
3943-c
3944-c Note nspin=1 in call to exchng...
3945-c
3946- CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
3947- FX = FX + DS(IS) * EXUNIF * F
3948-
3949- DKFDD = THD * KFS / DS(IS)
3950- DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
3951- DF1DD = 2 * (F1-1) * DSDD / S
3952- DFDD = KAPPA * DF1DD / F1**2
3953- DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
3954-
3955- DO 50 IX = 1,3
3956- GDS = 2 * GD(IX,IS)
3957- DSDGD = (S / GDMS) * GDS / GDMS
3958- DF1DGD = 2 * MU * S * DSDGD / KAPPA
3959- DFDGD = KAPPA * DF1DGD / F1**2
3960- DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
3961- 50 CONTINUE
3962- 60 CONTINUE
3963- FX = HALF * FX / DT
3964-
3965-C Set output arguments
3966- EX = FX
3967- EC = FC
3968- DO 90 IS = 1,nspin
3969- DEXDD(IS) = DFXDD(IS)
3970- DECDD(IS) = DFCDD(IS)
3971- DO 80 IX = 1,3
3972- DEXDGD(IX,IS) = DFXDGD(IX,IS)
3973- DECDGD(IX,IS) = DFCDGD(IX,IS)
3974- 80 CONTINUE
3975- 90 CONTINUE
3976-
3977- END
3978-
3979
3980=== modified file 'Util/Gen-basis/Makefile'
3981--- Util/Gen-basis/Makefile 2018-10-16 11:10:27 +0000
3982+++ Util/Gen-basis/Makefile 2019-01-30 13:43:54 +0000
3983@@ -53,7 +53,7 @@
3984 basis_specs.o atom.o memoryinfo.o memory.o periodic_table.o\
3985 pseudopotential.o pxf.o dot.o atom_options.o ldau_specs.o arw.o \
3986 timer.o xml.o m_walltime.o read_xc_info.o gen-basis.o $(SYSOBJ) \
3987- xc.o bsc_xcmod.o m_cite.o local_sys.o
3988+ m_cite.o local_sys.o
3989
3990 OBJS_IONCAT=f2kcli.o m_getopts.o basis_types.o precision.o parallel.o \
3991 parsing.o m_io.o io.o alloc.o atom_options.o basis_io.o atm_transfer.o atm_types.o \
3992@@ -61,7 +61,7 @@
3993 pseudopotential.o chemical.o m_filter.o \
3994 basis_specs.o atom.o memoryinfo.o m_memory.o periodic_table.o pxf.o \
3995 atmfuncs.o spher_harm.o interpolation.o old_atmfuncs.o arw.o \
3996- local_sys.o xc.o bsc_xcmod.o xml.o m_walltime.o ioncat.o ldau_specs.o m_cite.o $(SYSOBJ)
3997+ local_sys.o xml.o m_walltime.o ioncat.o ldau_specs.o m_cite.o $(SYSOBJ)
3998
3999
4000 #
4001@@ -138,9 +138,9 @@
4002 atm_transfer.o: periodic_table.o radial.o
4003 atm_types.o: precision.o radial.o
4004 atmfuncs.o: atm_types.o local_sys.o precision.o radial.o spher_harm.o
4005-atom.o: atmparams.o atom_options.o basis_specs.o basis_types.o bsc_xcmod.o
4006-atom.o: interpolation.o local_sys.o m_filter.o old_atmfuncs.o periodic_table.o
4007-atom.o: precision.o pseudopotential.o
4008+atom.o: atmparams.o atom_options.o basis_specs.o basis_types.o interpolation.o
4009+atom.o: local_sys.o m_filter.o old_atmfuncs.o periodic_table.o precision.o
4010+atom.o: pseudopotential.o
4011 atom_graph.o: alloc.o atm_types.o atmfuncs.o class_OrbitalDistribution.o
4012 atom_graph.o: class_SpData1D.o class_SpData2D.o class_SpData2D.o
4013 atom_graph.o: class_Sparsity.o intrinsic_missing.o ldau_specs.o local_sys.o
4014@@ -169,16 +169,14 @@
4015 broadcast_projections.o: m_trialorbitalclass.o parallel.o siesta2wannier90.o
4016 broyden_optim.o: local_sys.o m_broyddj_nocomm.o m_memory.o parallel.o
4017 broyden_optim.o: precision.o siesta_options.o units.o
4018-bsc_cellxc.o: alloc.o bsc_xcmod.o cellxc_mod.o local_sys.o mesh.o
4019-bsc_cellxc.o: moremeshsubs.o parallel.o parallelsubs.o precision.o
4020-bsc_xcmod.o: local_sys.o parallel.o precision.o
4021+bsc_cellxc.o: alloc.o local_sys.o mesh.o moremeshsubs.o parallel.o
4022+bsc_cellxc.o: parallelsubs.o precision.o
4023 cart2frac.o: local_sys.o precision.o
4024 cell_broyden_optim.o: local_sys.o m_broyddj_nocomm.o parallel.o precision.o
4025 cell_broyden_optim.o: units.o zmatrix.o
4026 cell_fire_optim.o: alloc.o local_sys.o m_fire.o m_memory.o parallel.o
4027 cell_fire_optim.o: precision.o siesta_options.o zmatrix.o
4028 cellsubs.o: precision.o
4029-cellxc_mod.o: bsc_xcmod.o local_sys.o
4030 cgvc.o: alloc.o conjgr_old.o local_sys.o m_mpi_utils.o parallel.o precision.o
4031 cgvc.o: units.o
4032 cgvc_zmatrix.o: alloc.o conjgr.o local_sys.o m_mpi_utils.o parallel.o
4033@@ -239,13 +237,13 @@
4034 detover.o: alloc.o linpack.o parallel.o parallelsubs.o precision.o
4035 dfscf.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o local_sys.o mesh.o
4036 dfscf.o: meshdscf.o meshphi.o parallel.o parallelsubs.o precision.o
4037-dhscf.o: alloc.o atmfuncs.o bsc_xcmod.o cellxc_mod.o delk.o dfscf.o
4038-dhscf.o: doping_uniform.o files.o forhar.o iogrid_netcdf.o local_sys.o
4039-dhscf.o: m_charge_add.o m_efield.o m_hartree_add.o m_iorho.o m_mesh_node.o
4040-dhscf.o: m_ncdf_io.o m_ncdf_siesta.o m_partial_charges.o m_rhog.o m_spin.o
4041-dhscf.o: m_ts_global_vars.o m_ts_hartree.o m_ts_options.o m_ts_voltage.o mesh.o
4042-dhscf.o: meshdscf.o meshsubs.o moremeshsubs.o parallel.o parsing.o precision.o
4043-dhscf.o: rhofft.o rhoofd.o siesta_options.o units.o vmat.o
4044+dhscf.o: alloc.o atmfuncs.o delk.o dfscf.o doping_uniform.o files.o forhar.o
4045+dhscf.o: iogrid_netcdf.o local_sys.o m_charge_add.o m_efield.o m_hartree_add.o
4046+dhscf.o: m_iorho.o m_mesh_node.o m_ncdf_io.o m_ncdf_siesta.o
4047+dhscf.o: m_partial_charges.o m_rhog.o m_spin.o m_ts_global_vars.o
4048+dhscf.o: m_ts_hartree.o m_ts_options.o m_ts_voltage.o mesh.o meshdscf.o
4049+dhscf.o: meshsubs.o moremeshsubs.o parallel.o parsing.o precision.o rhofft.o
4050+dhscf.o: rhoofd.o siesta_options.o units.o vmat.o
4051 diag.o: alloc.o diag_option.o local_sys.o parallel.o precision.o
4052 diag2g.o: fermid.o intrinsic_missing.o local_sys.o parallel.o parallelsubs.o
4053 diag2g.o: precision.o
4054@@ -295,13 +293,13 @@
4055 fft.o: alloc.o fft1d.o local_sys.o m_timer.o mesh.o parallel.o parallelsubs.o
4056 fft.o: precision.o
4057 fft1d.o: local_sys.o parallel.o precision.o
4058-final_H_f_stress.o: alloc.o atomlist.o class_SpData2D.o compute_max_diff.o
4059-final_H_f_stress.o: dnaefs.o files.o grdsam.o kinefsm.o ldau.o ldau_specs.o
4060-final_H_f_stress.o: local_sys.o m_dipol.o m_energies.o m_forces.o m_hsx.o
4061-final_H_f_stress.o: m_mpi_utils.o m_ncdf_siesta.o m_ntm.o m_spin.o m_stress.o
4062-final_H_f_stress.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o metaforce.o
4063-final_H_f_stress.o: molecularmechanics.o naefs.o nlefsm.o overfsm.o parallel.o
4064-final_H_f_stress.o: siesta_geom.o siesta_options.o sparse_matrices.o
4065+final_H_f_stress.o: alloc.o atomlist.o class_SpData1D.o class_SpData2D.o
4066+final_H_f_stress.o: compute_max_diff.o dnaefs.o files.o grdsam.o kinefsm.o
4067+final_H_f_stress.o: ldau.o ldau_specs.o local_sys.o m_dipol.o m_energies.o
4068+final_H_f_stress.o: m_forces.o m_hsx.o m_mpi_utils.o m_ncdf_siesta.o m_ntm.o
4069+final_H_f_stress.o: m_spin.o m_stress.o m_ts_io.o m_ts_kpoints.o m_ts_options.o
4070+final_H_f_stress.o: metaforce.o molecularmechanics.o naefs.o nlefsm.o overfsm.o
4071+final_H_f_stress.o: parallel.o siesta_geom.o siesta_options.o sparse_matrices.o
4072 final_H_f_stress.o: spinorbit.o units.o
4073 find_kgrid.o: alloc.o minvec.o parallel.o precision.o units.o
4074 fire_optim.o: alloc.o m_fire.o m_mpi_utils.o parallel.o precision.o
4075@@ -678,8 +676,8 @@
4076 meshdscf.o: alloc.o atomlist.o m_dscfcomm.o meshphi.o parallel.o parallelsubs.o
4077 meshdscf.o: precision.o
4078 meshphi.o: alloc.o precision.o
4079-meshsubs.o: alloc.o atm_types.o atmfuncs.o cellsubs.o cellxc_mod.o chkgmx.o
4080-meshsubs.o: fft1d.o local_sys.o mesh.o meshphi.o moremeshsubs.o parallel.o
4081+meshsubs.o: alloc.o atm_types.o atmfuncs.o cellsubs.o chkgmx.o fft1d.o
4082+meshsubs.o: local_sys.o mesh.o meshphi.o moremeshsubs.o parallel.o
4083 meshsubs.o: parallelsubs.o precision.o radial.o siesta_cml.o
4084 metaforce.o: alloc.o local_sys.o parallel.o precision.o
4085 minvec.o: cellsubs.o local_sys.o precision.o sorting.o
4086@@ -777,9 +775,9 @@
4087 rhoofd.o: precision.o
4088 rhoofdsp.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o local_sys.o
4089 rhoofdsp.o: mesh.o meshdscf.o meshphi.o precision.o
4090-save_density_matrix.o: atomlist.o files.o iodm_netcdf.o m_energies.o m_iodm.o
4091-save_density_matrix.o: m_matio.o m_ncdf_siesta.o m_spin.o m_steps.o
4092-save_density_matrix.o: m_ts_global_vars.o m_ts_iodm.o m_ts_options.o
4093+save_density_matrix.o: atomlist.o class_SpData2D.o files.o iodm_netcdf.o
4094+save_density_matrix.o: m_energies.o m_iodm.o m_matio.o m_ncdf_siesta.o m_spin.o
4095+save_density_matrix.o: m_steps.o m_ts_global_vars.o m_ts_iodm.o m_ts_options.o
4096 save_density_matrix.o: precision.o siesta_geom.o siesta_options.o
4097 save_density_matrix.o: sparse_matrices.o
4098 savepsi.o: alloc.o parallel.o parallelsubs.o precision.o
4099@@ -842,17 +840,16 @@
4100 siesta_forces.o: siesta_master.o siesta_options.o sparse_matrices.o
4101 siesta_forces.o: state_analysis.o state_init.o timer.o units.o write_subs.o
4102 siesta_geom.o: precision.o
4103-siesta_init.o: alloc.o atomlist.o bands.o bsc_xcmod.o
4104-siesta_init.o: class_Fstack_Pair_Geometry_SpData2D.o diag_option.o files.o
4105-siesta_init.o: flook_siesta.o ioxv.o kpoint_grid.o kpoint_pdos.o ksvinit.o
4106-siesta_init.o: local_sys.o m_check_walltime.o m_cite.o m_energies.o m_eo.o
4107-siesta_init.o: m_fixed.o m_forces.o m_iostruct.o m_mpi_utils.o m_new_dm.o
4108-siesta_init.o: m_rmaxh.o m_spin.o m_steps.o m_supercell.o m_timer.o
4109-siesta_init.o: m_wallclock.o metaforce.o molecularmechanics.o object_debug.o
4110-siesta_init.o: parallel.o parallelsubs.o projected_DOS.o siesta_cmlsubs.o
4111-siesta_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
4112-siesta_init.o: struct_init.o timer.o timestamp.o ts_init.o units.o writewave.o
4113-siesta_init.o: zmatrix.o
4114+siesta_init.o: alloc.o atomlist.o bands.o class_Fstack_Pair_Geometry_SpData2D.o
4115+siesta_init.o: diag_option.o files.o flook_siesta.o ioxv.o kpoint_grid.o
4116+siesta_init.o: kpoint_pdos.o ksvinit.o local_sys.o m_check_walltime.o m_cite.o
4117+siesta_init.o: m_energies.o m_eo.o m_fixed.o m_forces.o m_iostruct.o
4118+siesta_init.o: m_mpi_utils.o m_new_dm.o m_rmaxh.o m_spin.o m_steps.o
4119+siesta_init.o: m_supercell.o m_timer.o m_wallclock.o metaforce.o
4120+siesta_init.o: molecularmechanics.o object_debug.o parallel.o parallelsubs.o
4121+siesta_init.o: projected_DOS.o siesta_cmlsubs.o siesta_dicts.o siesta_geom.o
4122+siesta_init.o: siesta_options.o sparse_matrices.o struct_init.o timer.o
4123+siesta_init.o: timestamp.o ts_init.o units.o writewave.o zmatrix.o
4124 siesta_master.o: iopipes.o iosockets.o local_sys.o precision.o
4125 siesta_move.o: atomlist.o broyden_optim.o cell_broyden_optim.o
4126 siesta_move.o: cell_fire_optim.o dynamics.o fire_optim.o flook_siesta.o ioxv.o
4127@@ -926,7 +923,6 @@
4128 writewave.o: alloc.o atmfuncs.o atomlist.o densematrix.o diag.o diag_option.o
4129 writewave.o: get_kpoints_scale.o kpoint_grid.o local_sys.o m_spin.o parallel.o
4130 writewave.o: parallelsubs.o precision.o siesta_geom.o units.o
4131-xc.o: alloc.o bsc_xcmod.o local_sys.o precision.o
4132 xml.o: precision.o
4133 zm_broyden_optim.o: local_sys.o m_broyddj_nocomm.o parallel.o precision.o
4134 zm_broyden_optim.o: units.o zmatrix.o
4135
4136=== modified file 'version.info'
4137--- version.info 2019-01-30 11:09:09 +0000
4138+++ version.info 2019-01-30 13:43:54 +0000
4139@@ -1,1 +1,7 @@
4140+<<<<<<< TREE
4141 siesta-4.1-1059
4142+=======
4143+siesta-4.1-1050--xc-4
4144+
4145+
4146+>>>>>>> MERGE-SOURCE

Subscribers

People subscribed via source and target branches