Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to deal with vectors in 3-D real space.
10 : ! **************************************************************************************************
11 :
12 : MODULE negf_vectors
13 : USE kinds, ONLY: dp
14 : USE particle_types, ONLY: particle_type
15 : USE qs_subsys_types, ONLY: qs_subsys_get,&
16 : qs_subsys_type
17 : #include "./base/base_uses.f90"
18 :
19 : IMPLICIT NONE
20 : PRIVATE
21 :
22 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_vectors'
23 :
24 : PUBLIC :: contact_direction_vector, projection_on_direction_vector
25 :
26 : CONTAINS
27 :
28 : ! **************************************************************************************************
29 : !> \brief compute direction vector of the given contact
30 : !> \param origin origin
31 : !> \param direction_vector direction vector
32 : !> \param origin_bias origin which will be used to apply the external bias
33 : !> (in contrast with 'origin' it does not include screening region)
34 : !> \param direction_vector_bias direction vector which will be used to apply the external bias
35 : !> (together with 'origin_bias' it defines a contact region where
36 : !> the external potential is kept constant)
37 : !> \param atomlist_screening atoms belonging to the contact's screening region
38 : !> \param atomlist_bulk atoms belonging to the contact's bulk region
39 : !> \param subsys QuickStep subsystem
40 : !> \author Sergey Chulkov
41 : ! **************************************************************************************************
42 16 : SUBROUTINE contact_direction_vector(origin, direction_vector, origin_bias, direction_vector_bias, &
43 8 : atomlist_screening, atomlist_bulk, subsys)
44 : REAL(kind=dp), DIMENSION(3), INTENT(out) :: origin, direction_vector, origin_bias, &
45 : direction_vector_bias
46 : INTEGER, DIMENSION(:), INTENT(in) :: atomlist_screening, atomlist_bulk
47 : TYPE(qs_subsys_type), POINTER :: subsys
48 :
49 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contact_direction_vector'
50 :
51 : INTEGER :: handle, iatom, natoms_bulk, &
52 : natoms_screening, nparticles
53 : REAL(kind=dp) :: proj, proj_max, proj_min, proj_min_bias
54 : REAL(kind=dp), DIMENSION(3) :: vector
55 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
56 :
57 8 : CALL timeset(routineN, handle)
58 8 : CALL qs_subsys_get(subsys, particle_set=particle_set)
59 :
60 8 : natoms_screening = SIZE(atomlist_screening)
61 8 : natoms_bulk = SIZE(atomlist_bulk)
62 8 : nparticles = SIZE(particle_set)
63 :
64 8 : CPASSERT(natoms_screening > 0)
65 8 : CPASSERT(natoms_bulk > 0)
66 8 : CPASSERT(nparticles > 0)
67 :
68 : ! geometrical centre of the first contact unit cell
69 32 : origin = particle_set(atomlist_screening(1))%r
70 32 : DO iatom = 2, natoms_screening
71 104 : origin = origin + particle_set(atomlist_screening(iatom))%r
72 : END DO
73 32 : origin = origin/REAL(natoms_screening, kind=dp)
74 :
75 : ! geometrical centre of the second contact unit cell
76 32 : direction_vector = particle_set(atomlist_bulk(1))%r
77 64 : DO iatom = 2, natoms_bulk
78 232 : direction_vector = direction_vector + particle_set(atomlist_bulk(iatom))%r
79 : END DO
80 32 : direction_vector = direction_vector/REAL(natoms_bulk, kind=dp)
81 :
82 : ! vector between the geometrical centers of the first and the second contact unit cells
83 32 : direction_vector = direction_vector - origin
84 :
85 : ! the point 'center_of_coords0' belongs to the first unit cell, so the lowest projection of any point
86 : ! from the first unit cell on the direction vector 'center_of_coords1 - center_of_coords0' should be <= 0
87 : proj_min = 0.0_dp
88 : proj_min_bias = 0.0_dp
89 40 : DO iatom = 1, natoms_screening
90 128 : vector = particle_set(atomlist_screening(iatom))%r - origin
91 32 : proj = projection_on_direction_vector(vector, direction_vector)
92 :
93 32 : IF (proj < proj_min) proj_min = proj
94 40 : IF (proj > proj_min_bias) proj_min_bias = proj
95 : END DO
96 :
97 : ! the point 'center_of_coords1' belongs to the given contact, so the highest projection should be >= 1
98 : proj_max = 1.0_dp
99 232 : DO iatom = 1, nparticles
100 896 : vector = particle_set(iatom)%r - origin
101 224 : proj = projection_on_direction_vector(vector, direction_vector)
102 :
103 232 : IF (proj > proj_max) proj_max = proj
104 : END DO
105 :
106 : ! adjust the origin, so it lies on a plane between the given contact and the scattering region
107 32 : origin_bias = origin + proj_min_bias*direction_vector
108 32 : origin = origin + proj_min*direction_vector
109 :
110 : ! rescale the vector, so the last atom of the given contact and the point 'origin + direction_vector' lie on
111 : ! the same plane parallel to the 'origin' plane -- which separates the contact from the scattering region.
112 32 : direction_vector_bias = (proj_max - proj_min_bias)*direction_vector
113 32 : direction_vector = (proj_max - proj_min)*direction_vector
114 :
115 8 : CALL timestop(handle)
116 8 : END SUBROUTINE contact_direction_vector
117 :
118 : ! **************************************************************************************************
119 : !> \brief project the 'vector' onto the direction 'vector0'. Both vectors should have the same origin.
120 : !> \param vector vector to project
121 : !> \param vector0 direction vector
122 : !> \return projection (in terms of the direction vector's length)
123 : !> \author Sergey Chulkov
124 : !> \note \parblock
125 : !> & vector
126 : !> /
127 : !> /
128 : !> /
129 : !> *---|------> vector0
130 : !> l = proj * | vector0 |
131 : !>\endparblock
132 : ! **************************************************************************************************
133 1642410 : PURE FUNCTION projection_on_direction_vector(vector, vector0) RESULT(proj)
134 : REAL(kind=dp), DIMENSION(3), INTENT(in) :: vector, vector0
135 : REAL(kind=dp) :: proj
136 :
137 : REAL(kind=dp) :: len2, len2_v0, len2_v1
138 : REAL(kind=dp), DIMENSION(3) :: vector1
139 :
140 6569640 : vector1 = vector - vector0
141 6569640 : len2 = DOT_PRODUCT(vector, vector)
142 6569640 : len2_v0 = DOT_PRODUCT(vector0, vector0)
143 6569640 : len2_v1 = DOT_PRODUCT(vector1, vector1)
144 :
145 1642410 : proj = 0.5_dp*((len2 - len2_v1)/len2_v0 + 1.0_dp)
146 1642410 : END FUNCTION projection_on_direction_vector
147 : END MODULE negf_vectors
|
