1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module td_calc_oct_m
22  use iso_c_binding
23  use forces_oct_m
24  use geometry_oct_m
25  use global_oct_m
26  use grid_oct_m
27  use hamiltonian_elec_base_oct_m
28  use hamiltonian_elec_oct_m
29  use lasers_oct_m
30  use loct_math_oct_m
31  use mesh_function_oct_m
32  use messages_oct_m
33  use mpi_oct_m
34  use namespace_oct_m
35  use profiling_oct_m
36  use states_elec_calc_oct_m
37  use states_elec_oct_m
38  use states_elec_dim_oct_m
39
40  implicit none
41
42  private
43  public ::       &
44    td_calc_tacc, &
45    td_calc_tvel, &
46    td_calc_ionch
47
48contains
49
50! ---------------------------------------------------------
51!> Electronic acceleration (to calculate harmonic spectrum...)
52!! It is calculated as:
53!!
54!! \f[
55!! d2<x>/dt2 = d<p>/dt + i<[H,[V_nl,x]]> =
56!!           = i<[V_l,p]> + i<[V_nl,p]> - E(t)N + i<[H,[V_nl,x]]>
57!! \f]
58!! \warning This subroutine only works if ions are not
59!!          allowed to move
60! ---------------------------------------------------------
61subroutine td_calc_tacc(namespace, gr, geo, st, hm, acc, time)
62  type(namespace_t),        intent(in)  :: namespace
63  type(grid_t),             intent(in)  :: gr
64  type(geometry_t),         intent(in)  :: geo
65  type(states_elec_t),      intent(in)  :: st
66  type(hamiltonian_elec_t), intent(in)  :: hm
67  FLOAT,                    intent(in)  :: time
68  FLOAT,                    intent(out) :: acc(MAX_DIM)
69
70  FLOAT :: field(MAX_DIM), x(MAX_DIM)
71  CMPLX, allocatable :: zpsi(:, :), hzpsi(:,:), hhzpsi(:,:), xzpsi(:,:,:), vnl_xzpsi(:,:)
72  integer  :: j, k, ik, ist, idim
73
74#if defined(HAVE_MPI)
75  FLOAT   :: y(MAX_DIM)
76#endif
77
78  PUSH_SUB(td_calc_tacc)
79
80  ! The term i<[V_l,p]> + i<[V_nl,p]> may be considered as equal but opposite to the
81  ! force exerted by the electrons on the ions. COMMENT: This has to be thought about.
82  ! Maybe we are forgetting something....
83  call total_force_calculate(gr, geo, hm%ep, st, acc, hm%lda_u_level)
84
85  ! Adds the laser contribution : i<[V_laser, p]>
86  ! WARNING: this ignores the possibility of non-electric td external fields.
87  field = M_ZERO
88  do j = 1, hm%ep%no_lasers
89    call laser_electric_field(hm%ep%lasers(j), field(1:gr%sb%dim), time, CNST(0.001))
90    acc(1:gr%mesh%sb%dim) = acc(1:gr%mesh%sb%dim) - st%qtot*field(1:gr%mesh%sb%dim)
91  end do
92
93  if(.not. hm%ep%non_local) then
94    POP_SUB(td_calc_tacc)
95    return
96  end if
97
98  ! And now, i<[H,[V_nl,x]]>
99  x = M_ZERO
100  SAFE_ALLOCATE(zpsi(1:gr%mesh%np_part, 1:st%d%dim))
101  SAFE_ALLOCATE(hzpsi (1:gr%mesh%np_part, 1:st%d%dim))
102  SAFE_ALLOCATE(hhzpsi(1:3, 1:gr%mesh%np_part))
103
104  do ik = st%d%kpt%start, st%d%kpt%end
105    do ist = st%st_start, st%st_end
106
107      call states_elec_get_state(st, gr%mesh, ist, ik, zpsi)
108
109      call zhamiltonian_elec_apply_single(hm, namespace, gr%mesh, zpsi, hzpsi, ist, ik)
110
111      SAFE_ALLOCATE(xzpsi    (1:gr%mesh%np_part, 1:st%d%dim, 1:3))
112      SAFE_ALLOCATE(vnl_xzpsi(1:gr%mesh%np_part, 1:st%d%dim))
113      xzpsi = M_z0
114      do k = 1, gr%mesh%np
115        do j = 1, gr%mesh%sb%dim
116          xzpsi(k, 1:st%d%dim, j) = gr%mesh%x(k, j)*zpsi(k, 1:st%d%dim)
117        end do
118      end do
119
120      do j = 1, gr%mesh%sb%dim
121        call zhamiltonian_elec_apply_single(hm, namespace, gr%mesh, xzpsi(:, :, j), vnl_xzpsi, ist, ik, &
122          terms = TERM_NON_LOCAL_POTENTIAL)
123
124        do idim = 1, st%d%dim
125          x(j) = x(j) - 2*st%occ(ist, ik)*TOFLOAT(zmf_dotp(gr%mesh, hzpsi(1:gr%mesh%np, idim), vnl_xzpsi(:, idim)))
126        end do
127      end do
128
129      xzpsi = M_z0
130      do k = 1, gr%mesh%np
131        do j = 1, gr%mesh%sb%dim
132          xzpsi(k, 1:st%d%dim, j) = gr%mesh%x(k, j)*hzpsi(k, 1:st%d%dim)
133        end do
134      end do
135
136      do j = 1, gr%mesh%sb%dim
137        call zhamiltonian_elec_apply_single(hm, namespace, gr%mesh, xzpsi(:, :, j), vnl_xzpsi, ist, ik, &
138          terms = TERM_NON_LOCAL_POTENTIAL)
139
140        do idim = 1, st%d%dim
141          x(j) = x(j) + 2*st%occ(ist, ik)*TOFLOAT(zmf_dotp(gr%mesh, zpsi(:, idim), vnl_xzpsi(:, idim)))
142        end do
143      end do
144      SAFE_DEALLOCATE_A(xzpsi)
145      SAFE_DEALLOCATE_A(vnl_xzpsi)
146
147    end do
148  end do
149  SAFE_DEALLOCATE_A(hzpsi)
150  SAFE_DEALLOCATE_A(hhzpsi)
151
152#if defined(HAVE_MPI)
153  if(st%parallel_in_states) then
154    call MPI_Allreduce(x(1), y(1), gr%mesh%sb%dim, MPI_FLOAT, MPI_SUM, st%mpi_grp%comm, mpi_err)
155    x = y
156  end if
157#endif
158  acc = acc + x
159
160  POP_SUB(td_calc_tacc)
161end subroutine td_calc_tacc
162
163! ---------------------------------------------------------
164!> Electronic velocity (to calculate harmonic spectrum...)
165!! It is calculated as:
166!! \f[
167!! d<x>/dt = <p>
168!! \f]
169! ---------------------------------------------------------
170subroutine td_calc_tvel(gr, st, vel)
171  type(grid_t),        intent(in)  :: gr
172  type(states_elec_t), intent(in)  :: st
173  FLOAT,               intent(out) :: vel(MAX_DIM)
174
175  FLOAT, allocatable :: momentum(:,:,:)
176
177  PUSH_SUB(td_calc_tvel)
178
179  SAFE_ALLOCATE(momentum(1:gr%mesh%sb%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
180  call states_elec_calc_momentum(st, gr%der, momentum)
181
182  momentum(1:gr%mesh%sb%dim, st%st_start:st%st_end, 1) = &
183    sum(momentum(1:gr%mesh%sb%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end), 3)
184  momentum(1:gr%mesh%sb%dim, 1, 1) = &
185    sum(momentum(1:gr%mesh%sb%dim, st%st_start:st%st_end, 1), 2)
186  vel = momentum(1:gr%mesh%sb%dim, 1, 1)
187
188  SAFE_DEALLOCATE_A(momentum)
189  POP_SUB(td_calc_tvel)
190end subroutine td_calc_tvel
191
192! ---------------------------------------------------------
193!> Multiple ionization probabilities calculated form the KS orbital densities
194!! C. Ullrich, Journal of Molecular Structure: THEOCHEM 501, 315 (2000).
195! ---------------------------------------------------------
196subroutine td_calc_ionch(gr, st, ch, Nch)
197  type(grid_t),        intent(in)    :: gr
198  type(states_elec_t), intent(in)    :: st
199  integer,             intent(in)    :: Nch
200  FLOAT,               intent(out)   :: ch(0:Nch)
201
202  integer :: ik, ist, ii, jj, idim, Nid
203  FLOAT   :: prod, prod0
204  FLOAT, allocatable :: N(:), Nnot(:)
205  CMPLX,  allocatable :: zpsi(:)
206  !combinations
207  integer :: next
208  type(c_ptr) :: c
209  integer, allocatable :: idx0(:), idx(:), idxref(:)
210
211#if defined(HAVE_MPI)
212  FLOAT, allocatable :: Nbuf(:)
213#endif
214
215
216  PUSH_SUB(td_calc_ionch)
217
218  SAFE_ALLOCATE(   N(1: Nch))
219  SAFE_ALLOCATE(Nnot(1: Nch))
220  SAFE_ALLOCATE(zpsi(1:gr%mesh%np))
221
222  N(:)   = M_ZERO
223  Nnot(:)= M_ZERO
224
225
226  ii = 1
227  do ik = 1, st%d%nik
228    do ist = 1, st%nst
229      do idim = 1, st%d%dim
230
231        if (st%st_start <= ist .and. ist <= st%st_end .and. &
232              st%d%kpt%start <= ik .and. ik <= st%d%kpt%end) then
233          call states_elec_get_state(st, gr%mesh, idim, ist, ik, zpsi)
234          N(ii) = zmf_integrate(gr%mesh, zpsi(:) * conjg(zpsi(:)) )
235          Nnot(ii) = M_ONE - N(ii)
236        end if
237
238        ii = ii + 1
239
240      end do
241    end do
242  end do
243
244#if defined(HAVE_MPI)
245
246
247  if(st%parallel_in_states) then
248    SAFE_ALLOCATE(Nbuf(1: Nch))
249    Nbuf(:) = M_ZERO
250    call MPI_Allreduce(N(1), Nbuf(1), Nch, MPI_FLOAT, MPI_SUM, st%mpi_grp%comm, mpi_err)
251    N(:) = Nbuf(:)
252
253    Nbuf(:) = M_ZERO
254    call MPI_Allreduce(Nnot(1), Nbuf(1), Nch, MPI_FLOAT, MPI_SUM, st%mpi_grp%comm, mpi_err)
255    Nnot(:) = Nbuf(:)
256
257    SAFE_DEALLOCATE_A(Nbuf)
258  end if
259#endif
260
261! print * ,mpi_world%rank, "N    =", N(:)
262! print * ,mpi_world%rank, "Nnot =", Nnot(:)
263
264  ch(:) = M_ZERO
265
266  Nid = Nch - 1
267  SAFE_ALLOCATE(idx(0:Nid))
268  SAFE_ALLOCATE(idxref(0:Nid))
269  idxref = (/ (ii, ii = 0, Nid) /)
270  do ii = 0, Nch
271!     print *, "Nch= ", Nch, "ii", ii
272    call loct_combination_init(c, Nch, ii)
273    if (ii == 0 ) then
274      SAFE_ALLOCATE(idx0(0:0))
275    else
276      SAFE_ALLOCATE(idx0(0:ii-1))
277    end if
278!     print *,"size(idx0,1)=",size(idx0,1)
279    if(debug%info) then
280      call messages_write("P(")
281      call messages_write(ii)
282      call messages_write(") = 0")
283      call messages_new_line()
284    end if
285    ! Loop over combinations
286    do
287      prod  = M_ONE
288      prod0 = M_ONE
289      if(debug%info) then
290        call messages_write("P(")
291        call messages_write(ii)
292        call messages_write(") += ")
293      end if
294      call loct_get_combination(c, idx0)
295      idx(:) = idxref(:)
296      do jj = 0, ii-1
297        idx(idx0(jj))= -1
298        if(debug%info) then
299          call messages_write(" No(")
300          call messages_write(idx0(jj)+1)
301          call messages_write(") ")
302        end if
303        prod0 = prod0 * Nnot(idx0(jj)+1)
304      end do
305
306      do jj = 0, Nid
307        if(idx(jj)>=0) then
308          if(debug%info) then
309            call messages_write(" N(")
310            call messages_write(idx(jj)+1)
311            call messages_write(") ")
312          end if
313          prod = prod * N(idx(jj)+1)
314        end if
315      end do
316
317      if(debug%info) call messages_new_line()
318
319      ch(ii) = ch(ii) + prod*prod0
320
321      call loct_combination_next(c, next)
322      if( next /=  0) exit
323    end do
324    SAFE_DEALLOCATE_A(idx0)
325    call loct_combination_end(c)
326
327    if(debug%info) call messages_info()
328  end do
329
330  SAFE_DEALLOCATE_A(idx)
331  SAFE_DEALLOCATE_A(idxref)
332
333
334  SAFE_DEALLOCATE_A(N)
335  SAFE_DEALLOCATE_A(Nnot)
336  SAFE_DEALLOCATE_A(zpsi)
337
338  POP_SUB(td_calc_ionch)
339end subroutine td_calc_ionch
340
341
342end module td_calc_oct_m
343
344!! Local Variables:
345!! mode: f90
346!! coding: utf-8
347!! End:
348