Line data Source code
1 : !!****m* ABINIT/m_pred_steepdesc
2 : !! NAME
3 : !! m_pred_steepdesc
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, SE)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_pred_steepdesc
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_abimover
26 : use m_abihist
27 :
28 : use m_geometry, only : mkradim, mkrdim, xcart2xred, xred2xcart
29 : use m_predtk, only : fdtion
30 :
31 : implicit none
32 :
33 : private
34 : !!***
35 :
36 : public :: pred_steepdesc
37 : !!***
38 :
39 : contains
40 : !!***
41 :
42 : !!****f* ABINIT/pred_steepdesc
43 : !! NAME
44 : !! pred_steepdesc
45 : !!
46 : !! FUNCTION
47 : !! Ionmov predictor (21) Steepest Descent Algorithm
48 : !! The update of positions is given by the following equation:
49 : !!
50 : !! $$\Delta r_{n,i}=\lambda F_{n,i}$$
51 : !!
52 : !! r is the position of the 'n' ion along the 'i' direction
53 : !! F is the force of the 'n' ion along the 'i' direction.
54 : !!
55 : !! INPUTS
56 : !! ab_mover<type abimover>=Subset of dtset only related with
57 : !! | movement of ions and acell, contains:
58 : !! | dtion: Time step
59 : !! ! natom: Number of atoms
60 : !! | vis: viscosity
61 : !! | iatfix: Index of atoms and directions fixed
62 : !! | amass: Mass of ions
63 : !! icycle: Index of the internal cycle inside a time step (itime)
64 : !! itime: Index of time iteration
65 : !! zDEBUG : if true print some debugging information
66 : !!
67 : !! OUTPUT
68 : !!
69 : !! SIDE EFFECTS
70 : !! hist<type abihist>=Historical record of positions, forces, stresses, cell and energies.
71 : !!
72 : !! ncycle: Number of cycles of a particular time step
73 : !!
74 : !! NOTES
75 : !! * This routine is a predictor, it only produces new positions
76 : !! to be computed in the next iteration, this routine should
77 : !! produce not output at all
78 : !!
79 : !! SOURCE
80 :
81 0 : subroutine pred_steepdesc(ab_mover,forstr,hist,itime,zDEBUG,iexit)
82 :
83 : !Arguments ------------------------------------
84 : !scalars
85 : type(abimover),intent(in) :: ab_mover
86 : type(abihist),intent(inout),target :: hist
87 : type(abiforstr),intent(in) :: forstr
88 : integer,intent(in) :: itime,iexit
89 : logical,intent(in) :: zDEBUG
90 :
91 : !Local variables-------------------------------
92 : !scalars
93 : integer :: kk,jj,ihist_prev
94 : real(dp) :: em
95 : real(dp) :: f_cart
96 : real(dp) :: xc,str
97 : real(dp) :: xnow,lambda
98 : real(dp),save :: hh
99 : !arrays
100 : real(dp) :: acell(3),strten(6)
101 : real(dp) :: rprim(3,3),rprimd(3,3)
102 0 : real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom)
103 0 : real(dp) :: residual(3,ab_mover%natom)
104 0 : real(dp), ABI_CONTIGUOUS pointer :: fcart(:,:),vel(:,:)
105 :
106 : !***************************************************************************
107 : !Beginning of executable session
108 : !***************************************************************************
109 :
110 0 : if(iexit/=0)then
111 : return
112 : end if
113 :
114 : !write(std_out,*) '01'
115 : !##########################################################
116 : !### 01. Copy from the history to the variables
117 0 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
118 :
119 0 : do jj=1,3
120 0 : rprim(jj,1:3)=rprimd(jj,1:3)/acell(1:3)
121 : end do
122 :
123 0 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
124 0 : strten(:)=hist%strten(:,hist%ihist)
125 0 : fcart => hist%fcart(:,:,hist%ihist)
126 0 : vel => hist%vel(:,:,hist%ihist)
127 :
128 : !Fill the residual with forces (No preconditioning)
129 : !Or the preconditioned forces
130 0 : if (ab_mover%goprecon==0)then
131 0 : residual(:,:)=fcart(:,:)
132 : else
133 0 : residual(:,:)= forstr%fcart(:,:)
134 : end if
135 :
136 : !write(std_out,*) '01'
137 : !##########################################################
138 : !### 02. Get or compute de time step dtion
139 :
140 0 : if (ab_mover%dtion>0)then
141 0 : hh = ab_mover%dtion
142 : else
143 0 : hh=fdtion(ab_mover,itime,xcart,fcart,vel)
144 : end if
145 :
146 0 : lambda=hh
147 0 : write(std_out,*) 'Lambda',lambda
148 :
149 : !write(std_out,*) '02'
150 : !##########################################################
151 : !### 02. For all atoms and directions
152 0 : do kk=1,ab_mover%natom
153 : ! Normally this is the mass of the atom
154 : ! em=ab_mover%amass(kk)
155 : ! But the steepest algorithm is using unitary mass
156 : em=1
157 0 : do jj=1,3
158 :
159 : ! write(std_out,*) '03'
160 : ! ##########################################################
161 : ! ### 03. Filling other values from history (forces and vel)
162 0 : f_cart=residual(jj,kk)
163 0 : xc=xcart(jj,kk)
164 : ! This lambda is taken from the kinematical equation
165 : ! lambda=hh*hh/(2*em)
166 :
167 : ! write(std_out,*) '04'
168 : ! ##########################################################
169 : ! ### 04. Take first the atoms that are not allowed to move along
170 : ! ### this direction
171 : ! ### Warning : implemented in cartesian coordinates
172 0 : if (ab_mover%iatfix(jj,kk)==1) then
173 : ! Their positions will be the same as xcart
174 : xnow=xc
175 : else
176 :
177 : ! This is the main expresion (1)
178 0 : xnow=xc+lambda*f_cart
179 :
180 : end if !if(ab_mover%iatfix(jj,kk)==1)
181 :
182 : ! write(std_out,*) '05'
183 : ! ##########################################################
184 : ! ### 08. Update history
185 :
186 0 : xcart(jj,kk)=xnow
187 :
188 : ! write(std_out,*) '06'
189 : ! ##########################################################
190 : ! ### 09. End loops of atoms and directions
191 : end do ! jj=1,3
192 : end do ! kk=1,ab_mover%natom
193 :
194 0 : if (ab_mover%optcell/=0)then
195 :
196 0 : if (ab_mover%optcell==1)then
197 0 : do jj=1,3
198 0 : acell(jj)=acell(jj)+lambda*strten(jj)
199 : end do ! jj=1,3
200 0 : call mkrdim(acell,rprim,rprimd)
201 0 : elseif (ab_mover%optcell==2)then
202 0 : do kk=1,3
203 0 : do jj=1,3
204 0 : if (jj==1 .and. kk==1) str=strten(1)
205 0 : if (jj==2 .and. kk==2) str=strten(2)
206 0 : if (jj==3 .and. kk==3) str=strten(3)
207 0 : if (jj==1 .and. kk==2) str=strten(6)
208 0 : if (jj==1 .and. kk==3) str=strten(5)
209 0 : if (jj==2 .and. kk==1) str=strten(6)
210 0 : if (jj==3 .and. kk==1) str=strten(5)
211 0 : if (jj==2 .and. kk==3) str=strten(4)
212 0 : if (jj==3 .and. kk==2) str=strten(4)
213 0 : rprimd(jj,kk)=rprimd(jj,kk)+lambda*str
214 : end do ! jj=1,3
215 : end do ! kk=1,3
216 0 : call mkradim(acell,rprim,rprimd)
217 : end if
218 :
219 : end if
220 :
221 :
222 : !write(std_out,*) '08'
223 : !##########################################################
224 : !### 10. Filling history with the new values
225 :
226 : !Increase indices
227 0 : hist%ihist = abihist_findIndex(hist,+1)
228 :
229 : !Compute xred from xcart, and rprimd
230 0 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
231 :
232 0 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
233 0 : ihist_prev = abihist_findIndex(hist,-1)
234 0 : hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
235 :
236 0 : end subroutine pred_steepdesc
237 : !!***
238 :
239 : end module m_pred_steepdesc
240 : !!***
|
