Line data Source code
1 : !!****m* ABINIT/m_predtk
2 : !! NAME
3 : !! m_predtk
4 : !!
5 : !! FUNCTION
6 : !! Low-level procedures used by 45_geomoptim routines
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, SE)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_predtk
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_abimover
27 :
28 : implicit none
29 :
30 : private
31 : !!***
32 :
33 : public :: fdtion
34 : public :: prtxvf ! Print the values of xcart, vel, and fcart to unit iout.
35 : !!***
36 :
37 : contains
38 : !!***
39 :
40 : !!****f* ABINIT/fdtion
41 : !! NAME
42 : !! fdtion
43 : !!
44 : !! FUNCTION
45 : !! Compute the apropiated "dtion" from the present values
46 : !! of forces, velocity and viscosity
47 : !!
48 : !! INPUTS (in)
49 : !! hist<type abihist>=Historical record of positions, forces
50 : !! | acell, stresses, and energies,
51 : !! ab_mover<type abimover>=Subset of dtset only related with
52 : !! | movement of ions and acell, contains:
53 : !! | dtion: Time step
54 : !! ! natom: Number of atoms
55 : !! | vis: viscosity
56 : !! | iatfix: Index of atoms and directions fixed
57 : !! | amass: Mass of ions
58 : !! itime: Index of time iteration
59 : !! xcart(3,natom)= cartesian coordinates of atoms
60 : !! fcart(3,natom)= forces in cartesian coordinates
61 : !! vel(3,natom)= velocities
62 : !!
63 : !! OUTPUT (out)
64 : !! fdtion = time step computed
65 : !!
66 : !! SOURCE
67 :
68 0 : function fdtion(ab_mover,itime,xcart,fcart,vel)
69 :
70 : !Arguments ---------------------------------------------
71 : !scalars
72 : type(abimover),intent(in) :: ab_mover
73 : integer,intent(in) :: itime
74 : real(dp) :: fdtion
75 : !arrays
76 : real(dp) :: xcart(:,:),fcart(:,:),vel(:,:)
77 :
78 : !Local variables ------------------------------
79 : !scalars
80 : integer :: jj,kk
81 : real(dp) :: max,min,val
82 : real(dp) :: ff,xc,vv,em
83 :
84 : !************************************************************************
85 :
86 0 : max=0
87 0 : min=1e6
88 :
89 0 : do kk=1,ab_mover%natom
90 0 : em=ab_mover%amass(kk)
91 0 : do jj=1,3
92 0 : ff =fcart(jj,kk)
93 0 : xc =xcart(jj,kk)
94 0 : vv=vel(jj,kk)
95 :
96 0 : if (vv>1e-8) then
97 0 : val=abs(1.0_dp/vv)
98 0 : write(std_out,*) 'vel',kk,jj,val
99 0 : if (val>max) max=val
100 0 : if (val<min) min=val
101 : end if
102 :
103 0 : if (ff>1e-8) then
104 0 : val=sqrt(abs(2*em/ff))
105 0 : write(std_out,*) 'forces',kk,jj,val,em,ff
106 0 : if (val>max) max=val
107 0 : if (val<min) min=val
108 : end if
109 :
110 : end do
111 :
112 : end do
113 :
114 0 : write(std_out,*) "DTION max=",max
115 0 : write(std_out,*) "DTION min=",min
116 :
117 0 : if (itime==1)then
118 0 : fdtion=min/10
119 : else
120 0 : fdtion=min/10
121 : end if
122 :
123 0 : end function fdtion
124 : !!***
125 :
126 : !!****f* ABINIT/prtxvf
127 : !!
128 : !! NAME
129 : !! prtxvf
130 : !!
131 : !! FUNCTION
132 : !! Print the values of xcart, vel, and fcart to unit iout.
133 : !! Also compute and print max and rms forces.
134 :
135 : !! INPUTS
136 : !! fcart(3,natom)=forces (hartree/bohr)
137 : !! iatfix(3,natom)=1 for frozen or fixed atom along specified direction, else 0
138 : !! iout=unit number for printing
139 : !! natom=number of atoms in unit cell.
140 : !! prtvel=1 to print velocities, else do not print them
141 : !! vel(3,natom)=velocities
142 : !! xcart(3,natom)=cartesian coordinates (bohr)
143 : !!
144 : !! OUTPUT
145 : !! (only writing)
146 : !!
147 : !! SOURCE
148 :
149 405 : subroutine prtxvf(fcart,gred,iatfix,iout,natom,prtvel,vel,xcart,xred)
150 :
151 : !Arguments ------------------------------------
152 : !scalars
153 : integer,intent(in) :: iout,natom,prtvel
154 : !arrays
155 : integer,intent(in) :: iatfix(3,natom)
156 : real(dp),intent(in) :: fcart(3,natom),gred(3,natom)
157 : real(dp),intent(in) :: xcart(3,natom),xred(3,natom)
158 : real(dp),intent(in) :: vel(3,natom)
159 : !Local variables-------------------------------
160 : !scalars
161 : integer :: iatom,mu,unfixd
162 : real(dp) :: fmax,frms,val_max,val_rms
163 : character(len=500) :: msg
164 :
165 : ! *****************************************************************
166 :
167 405 : write(msg, '(a)' ) ' Cartesian coordinates (xcart) [bohr]'
168 405 : call wrtout(iout,msg,'COLL')
169 2799 : do iatom=1,natom
170 2394 : write(msg, '(1p,3e22.14)' )xcart(:,iatom)
171 2799 : call wrtout(iout,msg,'COLL')
172 : end do
173 :
174 405 : write(msg, '(a)' ) ' Reduced coordinates (xred)'
175 405 : call wrtout(iout,msg,'COLL')
176 2799 : do iatom=1,natom
177 2394 : write(msg, '(1p,3e22.14)' )xred(:,iatom)
178 2799 : call wrtout(iout,msg,'COLL')
179 : end do
180 :
181 : !Compute max |f| and rms f, EXCLUDING the components determined by iatfix
182 :
183 405 : fmax=0.0_dp
184 405 : frms=0.0_dp
185 405 : unfixd=0
186 2799 : do iatom=1,natom
187 9981 : do mu=1,3
188 9576 : if (iatfix(mu,iatom) /= 1) then
189 4014 : unfixd=unfixd+1
190 4014 : frms=frms+fcart(mu,iatom)**2
191 4014 : fmax=max(fmax,abs(fcart(mu,iatom)))
192 : end if
193 : end do
194 : end do
195 405 : if ( unfixd /= 0 ) frms=sqrt(frms/dble(unfixd))
196 :
197 : write(msg, '(a,1p,2e12.5,a)' ) &
198 405 : & ' Cartesian forces (fcart) [Ha/bohr]; max,rms=',fmax,frms,' (free atoms)'
199 405 : call wrtout(iout,msg,'COLL')
200 2799 : do iatom=1,natom
201 2394 : write(msg, '(1p,3e22.14)' )fcart(:,iatom)
202 2799 : call wrtout(iout,msg,'COLL')
203 : end do
204 :
205 405 : write(msg, '(a)' ) ' Gradient of E wrt nuclear positions in reduced coordinates (gred)'
206 405 : call wrtout(iout,msg,'COLL')
207 2799 : do iatom=1,natom
208 2394 : write(msg, '(1p,3e22.14)' )gred(:,iatom)
209 2799 : call wrtout(iout,msg,'COLL')
210 : end do
211 :
212 405 : if (prtvel == 1) then
213 :
214 : ! Compute max |v| and rms v,
215 : ! EXCLUDING the components determined by iatfix
216 0 : val_max=0.0_dp
217 0 : val_rms=0.0_dp
218 0 : unfixd=0
219 0 : do iatom=1,natom
220 0 : do mu=1,3
221 0 : if (iatfix(mu,iatom) /= 1) then
222 0 : unfixd=unfixd+1
223 0 : val_rms=val_rms+vel(mu,iatom)**2
224 0 : val_max=max(val_max,abs(vel(mu,iatom)**2))
225 : end if
226 : end do
227 : end do
228 0 : if ( unfixd /= 0 ) val_rms=sqrt(val_rms/dble(unfixd))
229 :
230 :
231 0 : write(msg, '(a,1p,2e12.5,a)' ) ' Cartesian velocities (vel) [bohr*Ha/hbar]; max,rms=',&
232 0 : & sqrt(val_max),val_rms,' (free atoms)'
233 0 : call wrtout(iout,msg,'COLL')
234 0 : do iatom=1,natom
235 0 : write(msg, '(1p,3e22.14)' ) vel(:,iatom)
236 0 : call wrtout(iout,msg,'COLL')
237 : end do
238 : end if
239 :
240 405 : end subroutine prtxvf
241 : !!***
242 :
243 : end module m_predtk
244 : !!***
|
