Line data Source code
1 : MODULE m_vdWfleur_grimme
2 : !
3 : ! implements Grimmes D2 and D3 based on routines from V. Caciuc (`19)
4 : !
5 :
6 : !
7 : ! Types for vdW-forces
8 : !
9 : PRIVATE
10 : TYPE atom_data
11 : INTEGER, DIMENSION(:), POINTER :: nr_atom_type
12 : INTEGER, DIMENSION(:), POINTER :: atomic_number
13 : REAL, DIMENSION(:,:),POINTER :: coord_bravais
14 : REAL, DIMENSION(:,:),POINTER :: coord_cart
15 : END TYPE atom_data
16 :
17 : real,allocatable:: force_stored(:,:)
18 : real :: energy_stored
19 : PUBLIC :: vdw_fleur_grimme
20 :
21 : CONTAINS
22 0 : SUBROUTINE vdW_fleur_grimme(input,atoms,sym,cell,e_vdW,f_vdW)
23 : USE m_constants,only:oUnit,tpi_const
24 : USE m_types,only: t_atoms,t_cell,t_sym,t_input
25 : USE DFT_D2, ONLY: driver_DFT_D2
26 : USE DFT_D3, ONLY: driver_DFT_D3
27 :
28 : IMPLICIT NONE
29 : TYPE(t_input),INTENT(IN) :: input
30 : TYPE(t_atoms),INTENT(IN) :: atoms
31 : TYPE(t_cell),INTENT(IN) :: cell
32 : TYPE(t_sym),INTENT(IN) :: sym
33 :
34 : REAL, INTENT (OUT) :: e_vdW,f_vdW(:,:)
35 :
36 : INTEGER :: NrAtomType,nsize,max_cyc,iop
37 : INTEGER :: NrAtoms,i,j,i1,i2,na,na1
38 : INTEGER,SAVE:: irep(3)=0 !will be determined at first call and reused later
39 : REAL :: start,finish,delta
40 : REAL :: test(3,8),brmin(3),brmax(3),force_i(3),f_rot(3)
41 : LOGICAL l_D2,l_in_au
42 : TYPE(atom_data) :: atom,atom_new
43 0 : REAL, ALLOCATABLE :: ener(:),force_vdW(:,:),force_max(:,:)
44 :
45 0 : if (allocated(force_stored)) THEN
46 0 : e_vdw=energy_stored
47 0 : f_vdw=force_stored
48 : return !vdW contribution already calculated
49 : endif
50 0 : max_cyc=merge(40,1,all(irep==0)) ! max. tries for bigger supercells
51 0 : ALLOCATE ( ener(max_cyc) )
52 :
53 0 : l_D2 = .false.
54 0 : l_in_au = .true. ! .false.
55 :
56 0 : NrAtomType=atoms%ntype
57 0 : ALLOCATE( atom%nr_atom_type(NrAtomType) ) ! neq(atoms%ntype)
58 0 : atom%nr_atom_type(:) = atoms%neq(:)
59 0 : NrAtoms = atoms%nat
60 0 : ALLOCATE( atom%atomic_number(NrAtoms) ) ! like zatom(natd) not (atoms%ntype)
61 0 : ALLOCATE( atom%coord_bravais(3,NrAtoms) ) ! taual(3,natd)
62 0 : ALLOCATE( atom%coord_cart(3,NrAtoms) ) ! pos(3,natd)
63 0 : ALLOCATE( force_vdW(3,NrAtoms),force_max(NrAtoms,max_cyc) )
64 :
65 0 : na = 0
66 0 : DO i1=1,atoms%ntype
67 0 : DO i2 = 1,atoms%neq(i1)
68 0 : na = na + 1
69 0 : atom%atomic_number(na) = NINT(atoms%zatom(i1))
70 0 : atom%coord_bravais(:,na)=matmul(cell%bmat,atoms%pos(:,na))/tpi_const
71 0 : atom%coord_cart(:,na)= atoms%pos(:,na)
72 : ENDDO
73 : ENDDO
74 :
75 :
76 :
77 0 : IF (max_cyc > 1) THEN ! determine a supercell size where the vdW energy converges
78 :
79 0 : test(:,1) = 0.0
80 0 : test(:,2) = cell%amat(1,:) ! tr!
81 0 : test(:,3) = cell%amat(2,:)
82 0 : test(:,4) = cell%amat(3,:)
83 0 : test(:,5) = test(:,2) + test(:,3)
84 0 : test(:,6) = test(:,3) + test(:,4)
85 0 : test(:,7) = test(:,4) + test(:,2)
86 0 : test(:,8) = test(:,5) + test(:,4)
87 :
88 0 : brmin(:) = minval(test(:,1:8),2)
89 0 : brmax(:) = maxval(test(:,1:8),2)
90 :
91 0 : irep(:) = 45.0 / (brmax(:)-brmin(:))
92 0 : IF (input%film) irep(3) = 0
93 :
94 : ENDIF
95 :
96 0 : cyc: DO nsize = 1, max_cyc
97 : ! generate supercell for vdW
98 0 : CALL gener_new_cell(atom,cell,irep,atom_new)
99 : !
100 : IF (l_D2) THEN
101 : ! using DFT-D2 method
102 : CALL driver_DFT_D2(atom%coord_cart,atom%atomic_number,atom_new%coord_cart,atom_new%atomic_number)
103 : ENDIF
104 : !
105 : ! note that in the driver_DFT_D3 atom%coord_cart and atom_new%coord_cart
106 : ! are modified from Angstrom to au [they have INTENT(INOUT)]
107 : ! using DFT-D3 method
108 0 : CALL driver_DFT_D3( atom%coord_cart,atom%atomic_number,atom_new%coord_cart,atom_new%atomic_number,l_in_au,ener(nsize),force_vdW)
109 :
110 0 : force_max(:,nsize) = sqrt(force_vdW(1,:)**2 + force_vdW(2,:)**2+ force_vdW(3,:)**2)
111 :
112 0 : IF (max_cyc == 1) THEN ! data taken from file
113 0 : delta = 0.0
114 0 : EXIT cyc
115 : ELSE
116 0 : IF (nsize > 1) THEN
117 0 : delta = ener(nsize)-ener(nsize-1)
118 0 : WRITE(oUnit,*) 'Delta = ',delta,' eV'
119 0 : DO i1=1,NrAtoms
120 0 : WRITE (oUnit,'(a15,i4,a3,3f15.9)') 'Delta vdW force',i1,' : ',force_max(i1,nsize)-force_max(i1,nsize-1)
121 : ENDDO
122 0 : IF (abs(delta) < input%vdw_tol) EXIT cyc
123 : ENDIF
124 : ENDIF
125 :
126 0 : irep(1:2) = irep(1:2) + 1
127 0 : IF (.not.input%film) irep(3) = irep(3) + 1
128 : ENDDO cyc
129 :
130 0 : IF (abs(delta) > input%vdW_tol) THEN
131 0 : WRITE (oUnit,*) 'vdW did not converge with cell size!'
132 0 : e_vdW = 0.0
133 : ELSE
134 0 : WRITE (oUnit,'(a13,3i5,a4,f12.3,a5)') 'vdW converged',irep(:)
135 0 : e_vdW = ener(nsize)/27.21138386
136 0 : WRITE ( oUnit,8060) e_vdW
137 0 : DO i1=1,NrAtoms
138 0 : WRITE (oUnit,'(a15,i4,a3,6f15.9)') 'vdW force atom ',i1,' : ',force_vdW(:,i1),atom%coord_cart(:,i1)
139 : ENDDO
140 : ENDIF
141 : 8060 FORMAT (/,/,' ----> vdW (D3) energy=',t40,f20.10,' htr')
142 :
143 0 : na = 0
144 0 : DO i1=1,NrAtomType
145 0 : f_rot=0.0
146 0 : DO i2 = 1,atoms%neq(i1) ! here symmetrization should be done
147 0 : na = na + 1
148 0 : iop = sym%ngopr(na) ! invtab(ngopr(na))
149 0 : force_i=matmul(cell%bmat,force_vdW(:,na))/tpi_const
150 0 : f_rot=f_rot+matmul(real(sym%mrot(:,:,iop)),force_i)
151 : ENDDO
152 0 : f_rot=f_rot/atoms%neq(i1)
153 0 : na1 = na - atoms%neq(i1) + 1
154 0 : force_i(:) = f_rot(:)
155 0 : IF (sym%invarind(na1) > 1) THEN
156 0 : DO i2 = 2, sym%invarind(na1)
157 0 : iop = sym%invarop(na1,i2)
158 0 : force_i=force_i+matmul(sym%mrot(:,:,iop),f_rot)
159 : ENDDO
160 0 : force_i(:) = force_i(:) / sym%invarind(na1)
161 : ENDIF
162 :
163 0 : f_vdW(:,i1)=matmul(cell%amat,force_i)
164 : ENDDO
165 0 : force_stored=f_vdW !Automatic allocate!
166 0 : energy_stored=e_vdw
167 : !
168 0 : END SUBROUTINE vdW_fleur_grimme
169 : !
170 :
171 0 : SUBROUTINE gener_new_cell(atom,cell,irep,atom_new)
172 :
173 : USE m_types,only: t_cell
174 :
175 : IMPLICIT NONE
176 :
177 : INTEGER, intent(IN) :: irep(3)
178 : TYPE(t_cell),intent(IN) :: cell
179 : TYPE(atom_data),intent(IN) :: atom
180 : TYPE(atom_data),intent(OUT) :: atom_new
181 :
182 :
183 : INTEGER :: irepeat1,irepeat2,irepeat3
184 : INTEGER :: i1,i2,j1,j2,j3
185 : INTEGER :: iatom,iline
186 : INTEGER :: NrNewAtoms
187 : !
188 0 : irepeat1=irep(1)
189 0 : irepeat2=irep(2)
190 0 : irepeat3=irep(3)
191 :
192 : !
193 0 : NrNewAtoms=(2*irepeat1+1)*(2*irepeat2+1)*(2*irepeat3+1)*sum(atom%nr_atom_type)
194 :
195 : !
196 0 : allocate(atom_new%atomic_number(1:NrNewAtoms))
197 0 : allocate(atom_new%coord_cart(3,NrNewAtoms))
198 0 : allocate(atom_new%coord_bravais(3,NrNewAtoms))
199 : !
200 : !
201 0 : iatom=0
202 0 : iline=0
203 0 : do i1=1,size(atom%nr_atom_type) ! how many different types of atoms
204 0 : do i2=1,atom%nr_atom_type(i1) ! how many atoms for each atom type
205 : ! iline is used to count the number of atoms in the initial unit cell
206 0 : iline=iline+1
207 0 : do j1=-irepeat1,irepeat1
208 0 : do j2=-irepeat2,irepeat2
209 0 : do j3=-irepeat3,irepeat3
210 0 : iatom=iatom+1
211 0 : atom_new%atomic_number(iatom)=atom%atomic_number(iline)
212 0 : atom_new%coord_bravais(:,iatom)=atom%coord_bravais(:,iline)+(/real(j1),real(j2),real(j3)/)
213 0 : atom_new%coord_cart(:,iatom)=matmul(cell%amat,atom_new%coord_bravais(:,iatom))
214 : enddo
215 : enddo
216 : enddo
217 : enddo
218 : enddo
219 : !
220 : !
221 0 : END SUBROUTINE gener_new_cell
222 :
223 0 : END MODULE m_vdWfleur_grimme
|
