72#include "./base/base_uses.f90"
77 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'negf_env_types'
78 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .true.
88 REAL(kind=
dp),
DIMENSION(3) :: direction_vector = -1.0_dp, origin = -1.0_dp
89 REAL(kind=
dp),
DIMENSION(3) :: direction_vector_bias = -1.0_dp, origin_bias = -1.0_dp
92 INTEGER :: direction_axis = -1
94 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atomlist_cell0
96 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atomlist_cell1
99 DIMENSION(:) :: atom_map_cell0, atom_map_cell1
101 REAL(kind=
dp) :: homo_energy = -1.0_dp
106 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: rho_00, rho_01
118 DIMENSION(:) :: contacts
132 INTEGER :: mixing_method = -1
139 TYPE negf_atom_map_contact_type
141 END TYPE negf_atom_map_contact_type
156 SUBROUTINE negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
162 INTEGER,
INTENT(in) :: log_unit
164 CHARACTER(len=*),
PARAMETER :: routinen =
'negf_env_create'
166 CHARACTER(len=default_string_length) :: contact_str, force_env_str, &
168 INTEGER :: handle, icontact, in_use, n_force_env, &
170 LOGICAL :: do_kpoints
172 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks_kp, matrix_s_kp
175 TYPE(negf_atom_map_contact_type),
ALLOCATABLE, &
176 DIMENSION(:) :: map_contact
181 CALL timeset(routinen, handle)
184 NULLIFY (sub_force_env)
185 CALL force_env_get(force_env, in_use=in_use, qs_env=qs_env, sub_force_env=sub_force_env)
187 IF (
ASSOCIATED(sub_force_env))
THEN
188 n_force_env =
SIZE(sub_force_env)
194 DO icontact = 1, n_force_env
195 CALL force_env_get(sub_force_env(icontact)%force_env, in_use=in_use)
201 cpabort(
"Quickstep is required for NEGF run.")
205 ncontacts =
SIZE(negf_control%contacts)
207 DO icontact = 1, ncontacts
208 IF (negf_control%contacts(icontact)%force_env_index > n_force_env)
THEN
209 WRITE (contact_str,
'(I11)') icontact
210 WRITE (force_env_str,
'(I11)') negf_control%contacts(icontact)%force_env_index
211 WRITE (n_force_env_str,
'(I11)') n_force_env
213 CALL cp_abort(__location__, &
214 "Contact number "//trim(adjustl(contact_str))//
" is linked with the FORCE_EVAL section number "// &
215 trim(adjustl(force_env_str))//
", however only "//trim(adjustl(n_force_env_str))// &
216 " FORCE_EVAL sections have been found. Note that FORCE_EVAL sections are enumerated from 0"// &
217 " and that the primary (0-th) section must contain all the atoms.")
224 CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, &
225 matrix_s_kp=matrix_s_kp, matrix_ks_kp=matrix_ks_kp, &
226 subsys=subsys, v_hartree_rspace=v_hartree_rspace)
229 cpabort(
"k-points are currently not supported for device FORCE_EVAL")
233 ALLOCATE (negf_env%contacts(ncontacts))
234 ALLOCATE (map_contact(ncontacts))
236 DO icontact = 1, ncontacts
237 IF (negf_control%contacts(icontact)%force_env_index > 0)
THEN
238 CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
239 CALL get_qs_env(qs_env_contact, subsys=subsys_contact)
241 CALL negf_env_contact_init_maps(contact_env=negf_env%contacts(icontact), &
242 contact_control=negf_control%contacts(icontact), &
243 atom_map=map_contact(icontact)%atom_map, &
244 eps_geometry=negf_control%eps_geometry, &
245 subsys_device=subsys, &
246 subsys_contact=subsys_contact)
248 IF (negf_env%contacts(icontact)%direction_axis == 0)
THEN
249 WRITE (contact_str,
'(I11)') icontact
250 WRITE (force_env_str,
'(I11)') negf_control%contacts(icontact)%force_env_index
251 CALL cp_abort(__location__, &
252 "One lattice vector of the contact unit cell (FORCE_EVAL section "// &
253 trim(adjustl(force_env_str))//
") must be parallel to the direction of the contact "// &
254 trim(adjustl(contact_str))//
".")
260 DO icontact = 1, ncontacts
261 IF (negf_control%contacts(icontact)%force_env_index > 0)
THEN
263 WRITE (log_unit,
'(/,T2,A,T70,I11)')
"NEGF| Construct the Kohn-Sham matrix for the contact ", icontact
265 CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
267 CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
269 CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
270 qs_env_contact=qs_env_contact, matrix_s_device=matrix_s_kp(1, 1)%matrix)
275 CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
278 DO icontact = 1, ncontacts
279 IF (negf_control%contacts(icontact)%force_env_index <= 0)
THEN
280 CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
281 contact_control=negf_control%contacts(icontact), &
282 sub_env=sub_env, qs_env=qs_env, &
283 eps_geometry=negf_control%eps_geometry)
288 CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
292 NULLIFY (negf_env%mixing_storage)
295 CALL get_qs_env(qs_env, dft_control=dft_control)
296 ALLOCATE (negf_env%mixing_storage)
298 negf_env%mixing_method, dft_control%qs_control%cutoff)
300 CALL timestop(handle)
313 SUBROUTINE negf_env_contact_init_maps(contact_env, contact_control, atom_map, &
314 eps_geometry, subsys_device, subsys_contact)
318 DIMENSION(:),
INTENT(inout) :: atom_map
319 REAL(kind=
dp),
INTENT(in) :: eps_geometry
322 CHARACTER(LEN=*),
PARAMETER :: routinen =
'negf_env_contact_init_maps'
324 INTEGER :: handle, natoms
326 CALL timeset(routinen, handle)
329 contact_env%direction_vector, &
330 contact_env%origin_bias, &
331 contact_env%direction_vector_bias, &
332 contact_control%atomlist_screening, &
333 contact_control%atomlist_bulk, &
336 contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys_contact, eps_geometry)
338 IF (contact_env%direction_axis /= 0)
THEN
339 natoms =
SIZE(contact_control%atomlist_bulk)
340 ALLOCATE (atom_map(natoms))
344 atom_list=contact_control%atomlist_bulk, &
345 subsys_device=subsys_device, &
346 subsys_contact=subsys_contact, &
347 eps_geometry=eps_geometry)
351 CALL list_atoms_in_bulk_primary_unit_cell(atomlist_cell0=contact_env%atomlist_cell0, &
352 atom_map_cell0=contact_env%atom_map_cell0, &
353 atomlist_bulk=contact_control%atomlist_bulk, &
355 origin=contact_env%origin, &
356 direction_vector=contact_env%direction_vector, &
357 direction_axis=contact_env%direction_axis, &
358 subsys_device=subsys_device)
361 CALL list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1=contact_env%atomlist_cell1, &
362 atom_map_cell1=contact_env%atom_map_cell1, &
363 atomlist_bulk=contact_control%atomlist_bulk, &
365 origin=contact_env%origin, &
366 direction_vector=contact_env%direction_vector, &
367 direction_axis=contact_env%direction_axis, &
368 subsys_device=subsys_device)
371 CALL timestop(handle)
372 END SUBROUTINE negf_env_contact_init_maps
382 SUBROUTINE negf_env_contact_init_matrices(contact_env, sub_env, qs_env_contact, matrix_s_device)
388 CHARACTER(LEN=*),
PARAMETER :: routinen =
'negf_env_contact_init_matrices'
390 INTEGER :: handle, iatom, ispin, nao, natoms, &
392 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_list0, atom_list1
393 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: index_to_cell, is_same_cell
394 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
395 LOGICAL :: do_kpoints
397 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
405 CALL timeset(routinen, handle)
408 dft_control=dft_control, &
409 do_kpoints=do_kpoints, &
411 matrix_ks_kp=matrix_ks_kp, &
412 matrix_s_kp=matrix_s_kp, &
416 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
418 CALL negf_homo_energy_estimate(contact_env%homo_energy, qs_env_contact)
420 natoms =
SIZE(contact_env%atomlist_cell0)
421 ALLOCATE (atom_list0(natoms))
423 atom_list0(iatom) = contact_env%atom_map_cell0(iatom)%iatom
427 IF (sum(abs(contact_env%atom_map_cell0(iatom)%cell(:))) > 0)
THEN
428 cpabort(
"NEGF K-points are not currently supported")
432 cpassert(
SIZE(contact_env%atomlist_cell1) == natoms)
433 ALLOCATE (atom_list1(natoms))
435 atom_list1(iatom) = contact_env%atom_map_cell1(iatom)%iatom
438 nspins = dft_control%nspins
439 nimages = dft_control%nimages
444 ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
445 cell_to_index(0, 0, 0) = 1
448 ALLOCATE (index_to_cell(3, nimages))
450 IF (.NOT. do_kpoints)
DEALLOCATE (cell_to_index)
454 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=sub_env%blacs_env)
457 ALLOCATE (contact_env%s_00, contact_env%s_01)
462 ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
463 ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
473 NULLIFY (matrix_s_ref)
475 CALL dbcsr_copy(matrix_s_ref, matrix_s_kp(1, 1)%matrix)
478 ALLOCATE (is_same_cell(natoms, natoms))
481 matrix_device=matrix_s_device, &
482 atom_list=contact_env%atomlist_cell0, &
483 atom_map=contact_env%atom_map_cell0, &
488 fm_cell1=contact_env%s_01, &
489 direction_axis=contact_env%direction_axis, &
490 matrix_kp=matrix_s_kp(1, :), &
491 index_to_cell=index_to_cell, &
492 atom_list0=atom_list0, atom_list1=atom_list1, &
493 subsys=subsys, mpi_comm_global=para_env, &
494 is_same_cell=is_same_cell, matrix_ref=matrix_s_ref)
499 fm_cell1=contact_env%h_01(ispin), &
500 direction_axis=contact_env%direction_axis, &
501 matrix_kp=matrix_ks_kp(ispin, :), &
502 index_to_cell=index_to_cell, &
503 atom_list0=atom_list0, atom_list1=atom_list1, &
504 subsys=subsys, mpi_comm_global=para_env, &
505 is_same_cell=is_same_cell)
508 fm_cell1=contact_env%rho_01(ispin), &
509 direction_axis=contact_env%direction_axis, &
510 matrix_kp=rho_ao_kp(ispin, :), &
511 index_to_cell=index_to_cell, &
512 atom_list0=atom_list0, atom_list1=atom_list1, &
513 subsys=subsys, mpi_comm_global=para_env, &
514 is_same_cell=is_same_cell)
517 DEALLOCATE (is_same_cell)
519 DEALLOCATE (index_to_cell)
520 DEALLOCATE (atom_list0, atom_list1)
522 CALL timestop(handle)
523 END SUBROUTINE negf_env_contact_init_matrices
534 SUBROUTINE negf_env_contact_init_matrices_gamma(contact_env, contact_control, sub_env, qs_env, eps_geometry)
539 REAL(kind=
dp),
INTENT(in) :: eps_geometry
541 CHARACTER(LEN=*),
PARAMETER :: routinen =
'negf_env_contact_init_matrices_gamma'
543 INTEGER :: handle, iatom, icell, ispin, nao_c, &
545 LOGICAL :: do_kpoints
546 REAL(kind=
dp),
DIMENSION(2) :: r2_origin_cell
547 REAL(kind=
dp),
DIMENSION(3) :: direction_vector, origin
549 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
556 CALL timeset(routinen, handle)
559 dft_control=dft_control, &
560 do_kpoints=do_kpoints, &
561 matrix_ks_kp=matrix_ks_kp, &
562 matrix_s_kp=matrix_s_kp, &
566 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
569 CALL cp_abort(__location__, &
570 "K-points in device region have not been implemented yet.")
573 nspins = dft_control%nspins
577 CALL cp_abort(__location__, &
578 "Primary and secondary bulk contact cells should be identical "// &
579 "in terms of the number of atoms of each kind, and their basis sets. "// &
580 "No single atom, however, can be shared between these two cells.")
583 contact_env%homo_energy = 0.0_dp
586 contact_env%direction_vector, &
587 contact_env%origin_bias, &
588 contact_env%direction_vector_bias, &
589 contact_control%atomlist_screening, &
590 contact_control%atomlist_bulk, &
593 contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys, eps_geometry)
598 origin = particle_set(contact_control%atomlist_screening(1))%r
599 DO iatom = 2,
SIZE(contact_control%atomlist_screening)
600 origin = origin + particle_set(contact_control%atomlist_screening(iatom))%r
602 origin = origin/real(
SIZE(contact_control%atomlist_screening), kind=
dp)
605 direction_vector = particle_set(contact_control%atomlist_cell(icell)%vector(1))%r
606 DO iatom = 2,
SIZE(contact_control%atomlist_cell(icell)%vector)
607 direction_vector = direction_vector + particle_set(contact_control%atomlist_cell(icell)%vector(iatom))%r
609 direction_vector = direction_vector/real(
SIZE(contact_control%atomlist_cell(icell)%vector), kind=
dp)
610 direction_vector = direction_vector - origin
611 r2_origin_cell(icell) = dot_product(direction_vector, direction_vector)
614 IF (sqrt(abs(r2_origin_cell(1) - r2_origin_cell(2))) < eps_geometry)
THEN
617 CALL cp_abort(__location__, &
618 "Primary and secondary bulk contact cells should not overlap ")
619 ELSE IF (r2_origin_cell(1) < r2_origin_cell(2))
THEN
620 ALLOCATE (contact_env%atomlist_cell0(
SIZE(contact_control%atomlist_cell(1)%vector)))
621 contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(1)%vector(:)
622 ALLOCATE (contact_env%atomlist_cell1(
SIZE(contact_control%atomlist_cell(2)%vector)))
623 contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(2)%vector(:)
625 ALLOCATE (contact_env%atomlist_cell0(
SIZE(contact_control%atomlist_cell(2)%vector)))
626 contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(2)%vector(:)
627 ALLOCATE (contact_env%atomlist_cell1(
SIZE(contact_control%atomlist_cell(1)%vector)))
628 contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(1)%vector(:)
632 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_c, ncol_global=nao_c, context=sub_env%blacs_env)
633 ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
634 ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
641 ALLOCATE (contact_env%s_00, contact_env%s_01)
648 fm=contact_env%h_00(ispin), &
649 atomlist_row=contact_env%atomlist_cell0, &
650 atomlist_col=contact_env%atomlist_cell0, &
651 subsys=subsys, mpi_comm_global=para_env, &
652 do_upper_diag=.true., do_lower=.true.)
654 fm=contact_env%h_01(ispin), &
655 atomlist_row=contact_env%atomlist_cell0, &
656 atomlist_col=contact_env%atomlist_cell1, &
657 subsys=subsys, mpi_comm_global=para_env, &
658 do_upper_diag=.true., do_lower=.true.)
661 fm=contact_env%rho_00(ispin), &
662 atomlist_row=contact_env%atomlist_cell0, &
663 atomlist_col=contact_env%atomlist_cell0, &
664 subsys=subsys, mpi_comm_global=para_env, &
665 do_upper_diag=.true., do_lower=.true.)
667 fm=contact_env%rho_01(ispin), &
668 atomlist_row=contact_env%atomlist_cell0, &
669 atomlist_col=contact_env%atomlist_cell1, &
670 subsys=subsys, mpi_comm_global=para_env, &
671 do_upper_diag=.true., do_lower=.true.)
675 fm=contact_env%s_00, &
676 atomlist_row=contact_env%atomlist_cell0, &
677 atomlist_col=contact_env%atomlist_cell0, &
678 subsys=subsys, mpi_comm_global=para_env, &
679 do_upper_diag=.true., do_lower=.true.)
681 fm=contact_env%s_01, &
682 atomlist_row=contact_env%atomlist_cell0, &
683 atomlist_col=contact_env%atomlist_cell1, &
684 subsys=subsys, mpi_comm_global=para_env, &
685 do_upper_diag=.true., do_lower=.true.)
687 CALL timestop(handle)
688 END SUBROUTINE negf_env_contact_init_matrices_gamma
699 SUBROUTINE negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
705 CHARACTER(LEN=*),
PARAMETER :: routinen =
'negf_env_device_init_matrices'
707 INTEGER :: handle, icontact, ispin, nao_c, nao_s, &
709 LOGICAL :: do_kpoints
712 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks_kp, matrix_s_kp
720 CALL timeset(routinen, handle)
722 IF (
ALLOCATED(negf_control%atomlist_S_screening))
THEN
724 dft_control=dft_control, &
725 do_kpoints=do_kpoints, &
726 matrix_ks_kp=matrix_ks_kp, &
727 matrix_s_kp=matrix_s_kp, &
731 CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
734 CALL cp_abort(__location__, &
735 "K-points in device region have not been implemented yet.")
738 ncontacts =
SIZE(negf_control%contacts)
739 nspins = dft_control%nspins
745 NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct)
746 ALLOCATE (negf_env%h_s(nspins))
748 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env)
749 ALLOCATE (negf_env%s_s)
754 ALLOCATE (negf_env%v_hartree_s)
759 ALLOCATE (negf_env%h_sc(nspins, ncontacts), negf_env%s_sc(ncontacts))
760 DO icontact = 1, ncontacts
762 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_c, context=sub_env%blacs_env)
767 CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
776 fm=negf_env%h_s(ispin), &
777 atomlist_row=negf_control%atomlist_S_screening, &
778 atomlist_col=negf_control%atomlist_S_screening, &
779 subsys=subsys, mpi_comm_global=para_env, &
780 do_upper_diag=.true., do_lower=.true.)
785 atomlist_row=negf_control%atomlist_S_screening, &
786 atomlist_col=negf_control%atomlist_S_screening, &
787 subsys=subsys, mpi_comm_global=para_env, &
788 do_upper_diag=.true., do_lower=.true.)
791 NULLIFY (hmat%matrix)
793 CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix)
796 CALL pw_pool%create_pw(v_hartree)
797 CALL negf_env_init_v_hartree(v_hartree, negf_env%contacts, negf_control%contacts)
799 CALL integrate_v_rspace(v_rspace=v_hartree, hmat=hmat, qs_env=qs_env, &
800 calculate_forces=.false., compute_tau=.false., gapw=.false.)
802 CALL pw_pool%give_back_pw(v_hartree)
805 fm=negf_env%v_hartree_s, &
806 atomlist_row=negf_control%atomlist_S_screening, &
807 atomlist_col=negf_control%atomlist_S_screening, &
808 subsys=subsys, mpi_comm_global=para_env, &
809 do_upper_diag=.true., do_lower=.true.)
814 DO icontact = 1, ncontacts
817 fm=negf_env%h_sc(ispin, icontact), &
818 atomlist_row=negf_control%atomlist_S_screening, &
819 atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
820 subsys=subsys, mpi_comm_global=para_env, &
821 do_upper_diag=.true., do_lower=.true.)
825 fm=negf_env%s_sc(icontact), &
826 atomlist_row=negf_control%atomlist_S_screening, &
827 atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
828 subsys=subsys, mpi_comm_global=para_env, &
829 do_upper_diag=.true., do_lower=.true.)
833 CALL timestop(handle)
834 END SUBROUTINE negf_env_device_init_matrices
843 SUBROUTINE negf_env_init_v_hartree(v_hartree, contact_env, contact_control)
846 INTENT(in) :: contact_env
848 INTENT(in) :: contact_control
850 CHARACTER(len=*),
PARAMETER :: routinen =
'negf_env_init_v_hartree'
851 REAL(kind=
dp),
PARAMETER :: threshold = 16.0_dp*epsilon(0.0_dp)
853 INTEGER :: dx, dy, dz, handle, icontact, ix, iy, &
854 iz, lx, ly, lz, ncontacts, ux, uy, uz
855 REAL(kind=
dp) :: dvol, pot, proj, v1, v2
856 REAL(kind=
dp),
DIMENSION(3) :: dirvector_bias, point_coord, &
857 point_indices, vector
859 CALL timeset(routinen, handle)
861 ncontacts =
SIZE(contact_env)
862 cpassert(
SIZE(contact_control) == ncontacts)
863 cpassert(ncontacts == 2)
865 dirvector_bias = contact_env(2)%origin_bias - contact_env(1)%origin_bias
866 v1 = contact_control(1)%v_external
867 v2 = contact_control(2)%v_external
869 lx = v_hartree%pw_grid%bounds_local(1, 1)
870 ux = v_hartree%pw_grid%bounds_local(2, 1)
871 ly = v_hartree%pw_grid%bounds_local(1, 2)
872 uy = v_hartree%pw_grid%bounds_local(2, 2)
873 lz = v_hartree%pw_grid%bounds_local(1, 3)
874 uz = v_hartree%pw_grid%bounds_local(2, 3)
876 dx = v_hartree%pw_grid%npts(1)/2
877 dy = v_hartree%pw_grid%npts(2)/2
878 dz = v_hartree%pw_grid%npts(3)/2
880 dvol = v_hartree%pw_grid%dvol
883 point_indices(3) = real(iz + dz, kind=
dp)
885 point_indices(2) = real(iy + dy, kind=
dp)
888 point_indices(1) = real(ix + dx, kind=
dp)
889 point_coord(:) = matmul(v_hartree%pw_grid%dh, point_indices)
891 vector = point_coord - contact_env(1)%origin_bias
893 IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp)
THEN
897 IF (proj < 0.0_dp)
THEN
899 ELSE IF (proj > 1.0_dp)
THEN
902 pot = v1 + (v2 - v1)*proj
905 DO icontact = 1, ncontacts
906 vector = point_coord - contact_env(icontact)%origin_bias
909 IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp)
THEN
910 pot = contact_control(icontact)%v_external
916 v_hartree%array(ix, iy, iz) = pot*dvol
921 CALL timestop(handle)
922 END SUBROUTINE negf_env_init_v_hartree
933 FUNCTION contact_direction_axis(direction_vector, subsys_contact, eps_geometry)
RESULT(direction_axis)
934 REAL(kind=
dp),
DIMENSION(3),
INTENT(in) :: direction_vector
936 REAL(kind=
dp),
INTENT(in) :: eps_geometry
937 INTEGER :: direction_axis
940 REAL(kind=
dp),
DIMENSION(3) :: scaled
950 IF (abs(scaled(i)) > eps_geometry)
THEN
951 IF (scaled(i) > 0.0_dp)
THEN
961 IF (naxes /= 1) direction_axis = 0
962 END FUNCTION contact_direction_axis
971 SUBROUTINE negf_homo_energy_estimate(homo_energy, qs_env)
972 REAL(kind=
dp),
INTENT(out) :: homo_energy
975 CHARACTER(LEN=*),
PARAMETER :: routinen =
'negf_homo_energy_estimate'
976 INTEGER,
PARAMETER :: gamma_point = 1
978 INTEGER :: handle, homo, ikpgr, ikpoint, imo, &
979 ispin, kplocal, nmo, nspins
980 INTEGER,
DIMENSION(2) :: kp_range
981 LOGICAL :: do_kpoints
982 REAL(kind=
dp) :: my_homo_energy
983 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigenvalues
987 TYPE(
mo_set_type),
DIMENSION(:, :),
POINTER :: mos_kp
990 CALL timeset(routinen, handle)
991 my_homo_energy = 0.0_dp
993 CALL get_qs_env(qs_env, para_env=para_env, mos=mos, kpoints=kpoints, do_kpoints=do_kpoints)
996 CALL get_kpoint_info(kpoints, kp_env=kp_env, kp_range=kp_range, para_env_kp=para_env_kp)
999 IF (para_env_kp%mepos == 0 .AND. kp_range(1) <= gamma_point .AND. kp_range(2) >= gamma_point)
THEN
1000 kplocal = kp_range(2) - kp_range(1) + 1
1002 DO ikpgr = 1, kplocal
1003 CALL get_kpoint_env(kp_env(ikpgr)%kpoint_env, nkpoint=ikpoint, mos=mos_kp)
1005 IF (ikpoint == gamma_point)
THEN
1007 CALL get_mo_set(mos_kp(1, 1), homo=homo, eigenvalues=eigenvalues)
1009 my_homo_energy = eigenvalues(homo)
1015 CALL para_env%sum(my_homo_energy)
1022 cpwarn(
"It is strongly advised to use k-points within all contact FORCE_EVAL-s")
1026 spin_loop:
DO ispin = 1, nspins
1027 CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
1030 IF (eigenvalues(imo) /= 0.0_dp)
EXIT spin_loop
1035 cpabort(
"Orbital transformation (OT) for contact FORCE_EVAL-s is not supported")
1038 my_homo_energy = eigenvalues(homo)
1041 homo_energy = my_homo_energy
1042 CALL timestop(handle)
1043 END SUBROUTINE negf_homo_energy_estimate
1059 SUBROUTINE list_atoms_in_bulk_primary_unit_cell(atomlist_cell0, atom_map_cell0, atomlist_bulk, atom_map, &
1060 origin, direction_vector, direction_axis, subsys_device)
1061 INTEGER,
ALLOCATABLE,
DIMENSION(:),
INTENT(inout) :: atomlist_cell0
1063 DIMENSION(:),
INTENT(inout) :: atom_map_cell0
1064 INTEGER,
DIMENSION(:),
INTENT(in) :: atomlist_bulk
1066 REAL(kind=
dp),
DIMENSION(3),
INTENT(in) :: origin, direction_vector
1067 INTEGER,
INTENT(in) :: direction_axis
1070 CHARACTER(LEN=*),
PARAMETER :: routinen =
'list_atoms_in_bulk_primary_unit_cell'
1072 INTEGER :: atom_min, dir_axis_min, &
1073 direction_axis_abs, handle, iatom, &
1074 natoms_bulk, natoms_cell0
1075 REAL(kind=
dp) :: proj, proj_min
1076 REAL(kind=
dp),
DIMENSION(3) :: vector
1079 CALL timeset(routinen, handle)
1080 CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1082 natoms_bulk =
SIZE(atomlist_bulk)
1083 cpassert(
SIZE(atom_map, 1) == natoms_bulk)
1084 direction_axis_abs = abs(direction_axis)
1089 DO iatom = 1, natoms_bulk
1090 vector = particle_set(atomlist_bulk(iatom))%r - origin
1093 IF (proj < proj_min)
THEN
1099 dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1102 DO iatom = 1, natoms_bulk
1103 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) &
1104 natoms_cell0 = natoms_cell0 + 1
1107 ALLOCATE (atomlist_cell0(natoms_cell0))
1108 ALLOCATE (atom_map_cell0(natoms_cell0))
1111 DO iatom = 1, natoms_bulk
1112 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min)
THEN
1113 natoms_cell0 = natoms_cell0 + 1
1114 atomlist_cell0(natoms_cell0) = atomlist_bulk(iatom)
1115 atom_map_cell0(natoms_cell0) = atom_map(iatom)
1119 CALL timestop(handle)
1120 END SUBROUTINE list_atoms_in_bulk_primary_unit_cell
1139 SUBROUTINE list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1, atom_map_cell1, atomlist_bulk, atom_map, &
1140 origin, direction_vector, direction_axis, subsys_device)
1141 INTEGER,
ALLOCATABLE,
DIMENSION(:),
INTENT(inout) :: atomlist_cell1
1143 DIMENSION(:),
INTENT(inout) :: atom_map_cell1
1144 INTEGER,
DIMENSION(:),
INTENT(in) :: atomlist_bulk
1146 REAL(kind=
dp),
DIMENSION(3),
INTENT(in) :: origin, direction_vector
1147 INTEGER,
INTENT(in) :: direction_axis
1150 CHARACTER(LEN=*),
PARAMETER :: routinen =
'list_atoms_in_bulk_secondary_unit_cell'
1152 INTEGER :: atom_min, dir_axis_min, &
1153 direction_axis_abs, handle, iatom, &
1154 natoms_bulk, natoms_cell1, offset
1155 REAL(kind=
dp) :: proj, proj_min
1156 REAL(kind=
dp),
DIMENSION(3) :: vector
1159 CALL timeset(routinen, handle)
1160 CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1162 natoms_bulk =
SIZE(atomlist_bulk)
1163 cpassert(
SIZE(atom_map, 1) == natoms_bulk)
1164 direction_axis_abs = abs(direction_axis)
1165 offset = sign(1, direction_axis)
1170 DO iatom = 1, natoms_bulk
1171 vector = particle_set(atomlist_bulk(iatom))%r - origin
1174 IF (proj < proj_min)
THEN
1180 dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1183 DO iatom = 1, natoms_bulk
1184 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) &
1185 natoms_cell1 = natoms_cell1 + 1
1188 ALLOCATE (atomlist_cell1(natoms_cell1))
1189 ALLOCATE (atom_map_cell1(natoms_cell1))
1192 DO iatom = 1, natoms_bulk
1193 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset)
THEN
1194 natoms_cell1 = natoms_cell1 + 1
1195 atomlist_cell1(natoms_cell1) = atomlist_bulk(iatom)
1196 atom_map_cell1(natoms_cell1) = atom_map(iatom)
1197 atom_map_cell1(natoms_cell1)%cell(direction_axis_abs) = dir_axis_min
1201 CALL timestop(handle)
1202 END SUBROUTINE list_atoms_in_bulk_secondary_unit_cell
1213 CHARACTER(len=*),
PARAMETER :: routinen =
'negf_env_release'
1215 INTEGER :: handle, icontact
1217 CALL timeset(routinen, handle)
1219 IF (
ALLOCATED(negf_env%contacts))
THEN
1220 DO icontact =
SIZE(negf_env%contacts), 1, -1
1221 CALL negf_env_contact_release(negf_env%contacts(icontact))
1224 DEALLOCATE (negf_env%contacts)
1234 IF (
ASSOCIATED(negf_env%s_s))
THEN
1236 DEALLOCATE (negf_env%s_s)
1237 NULLIFY (negf_env%s_s)
1244 IF (
ASSOCIATED(negf_env%v_hartree_s))
THEN
1246 DEALLOCATE (negf_env%v_hartree_s)
1247 NULLIFY (negf_env%v_hartree_s)
1251 IF (
ASSOCIATED(negf_env%mixing_storage))
THEN
1253 DEALLOCATE (negf_env%mixing_storage)
1256 CALL timestop(handle)
1263 SUBROUTINE negf_env_contact_release(contact_env)
1266 CHARACTER(len=*),
PARAMETER :: routinen =
'negf_env_contact_release'
1270 CALL timeset(routinen, handle)
1285 IF (
ASSOCIATED(contact_env%s_00))
THEN
1287 DEALLOCATE (contact_env%s_00)
1288 NULLIFY (contact_env%s_00)
1292 IF (
ASSOCIATED(contact_env%s_01))
THEN
1294 DEALLOCATE (contact_env%s_01)
1295 NULLIFY (contact_env%s_01)
1298 IF (
ALLOCATED(contact_env%atomlist_cell0))
DEALLOCATE (contact_env%atomlist_cell0)
1299 IF (
ALLOCATED(contact_env%atomlist_cell1))
DEALLOCATE (contact_env%atomlist_cell1)
1300 IF (
ALLOCATED(contact_env%atom_map_cell0))
DEALLOCATE (contact_env%atom_map_cell0)
1302 CALL timestop(handle)
1303 END SUBROUTINE negf_env_contact_release
Handles all functions related to the CELL.
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_set(matrix, alpha)
...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
integer, parameter, public use_qs_force
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_env(kpoint_env, nkpoint, wkp, xkp, is_local, mos)
Get information from a single kpoint environment.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Interface to the message passing library MPI.
Map atoms between various force environments.
subroutine, public negf_map_atomic_indices(atom_map, atom_list, subsys_device, subsys_contact, eps_geometry)
Map atoms in the cell 'subsys_device' listed in 'atom_list' with the corresponding atoms in the cell ...
Input control types for NEGF based quantum transport calculations.
Environment for NEGF based quantum transport calculations.
subroutine, public negf_env_release(negf_env)
Release a NEGF environment variable.
subroutine, public negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
Helper routines to manipulate with matrices.
subroutine, public negf_reference_contact_matrix(matrix_contact, matrix_device, atom_list, atom_map, para_env)
Extract part of the DBCSR matrix based on selected atoms and copy it into another DBCSR matrix.
integer function, public number_of_atomic_orbitals(subsys, atom_list)
Compute the number of atomic orbitals of the given set of atoms.
subroutine, public invert_cell_to_index(cell_to_index, nimages, index_to_cell)
Invert cell_to_index mapping between unit cells and DBCSR matrix images.
subroutine, public negf_copy_contact_matrix(fm_cell0, fm_cell1, direction_axis, matrix_kp, index_to_cell, atom_list0, atom_list1, subsys, mpi_comm_global, is_same_cell, matrix_ref)
Driver routine to extract diagonal and off-diagonal blocks from a symmetric DBCSR matrix.
subroutine, public negf_copy_sym_dbcsr_to_fm_submat(matrix, fm, atomlist_row, atomlist_col, subsys, mpi_comm_global, do_upper_diag, do_lower)
Extract part of the DBCSR matrix based on selected atoms and copy it into a dense matrix.
Environment for NEGF based quantum transport calculations.
Routines to deal with vectors in 3-D real space.
pure real(kind=dp) function, public projection_on_direction_vector(vector, vector0)
project the 'vector' onto the direction 'vector0'. Both vectors should have the same origin.
subroutine, public contact_direction_vector(origin, direction_vector, origin_bias, direction_vector_bias, atomlist_screening, atomlist_bulk, subsys)
compute direction vector of the given contact
Define the data structure for the particle information.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
module that contains the definitions of the scf types
subroutine, public mixing_storage_release(mixing_store)
releases a mixing_storage
subroutine, public mixing_storage_create(mixing_store, mixing_section, mixing_method, ecut)
creates a mixing_storage
Utility subroutine for qs energy calculation.
subroutine, public qs_energies_init(qs_env, calc_forces)
Refactoring of qs_energies_scf. Driver routine for the initial setup and calculations for a qs energy...
Perform a QUICKSTEP wavefunction optimization (single point)
subroutine, public qs_energies(qs_env, consistent_energies, calc_forces)
Driver routine for QUICKSTEP single point wavefunction optimization.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Integrate single or product functions over a potential on a RS grid.
Definition and initialisation of the mo data type.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Type defining parameters related to the simulation cell.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
allows for the creation of an array of force_env
wrapper to abstract the force evaluation of the various methods
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Structure that maps the given atom in the sourse FORCE_EVAL section with another atom from the target...
Input parameters related to the NEGF run.
Parallel (sub)group environment.
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
keeps the density in various representations, keeping track of which ones are valid.