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