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