Branch data 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), 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 : : !!***
|
