Dear Bartek,
I am very pleased to get your e-mail and thanks a lot for your kind help!
Ok, I will try to modify the velocity-version and force-version according
to your suggestions!
But, I have still another question on this problem: Can the applied
external force or velocity affect the convergence of the SCF calculations? Any
experience?
Thanks again!
Best Regards!
Sincerely Yours
S.F.Li
-----Original Message-----
From: Bartek Szyja on behalf of Bartek Szyja
Sent: Fri 1/21/2011 7:12 AM
To: [email protected]
Subject: Re: [SIESTA-L] For help: How to pply a constant initial velocity or
external force on some atoms of a system.
On Fri, 2011-01-21 at 12:36 +0100, Riccardo Rurali wrote:
> > 2)Alternatively, can we apply external force on the atoms?
> > If it is okay, how to construct the input file? Or should we modify the
> > source code?
>
> You should modify the source code, but the developers have prepared a
> simple way for you to do it.
> Edit the file constr.f in the Src folder and add the kind of constraint
> that you like. If you add the line
>
> fa(3,132) = 0.5
>
> for instance, at each MD step the force along z of the atoms number 132
> will be set to 0.5.
>
> Don't forget to recompile Siesta and to specify the constraint:
>
> %block GeometryConstraints
> routine constr
> %endblock GeometryConstraints
The only problem with it is that the recompilation of the code is
necessary after each change you make to the external force. I have been
studying similar subject before and dealt with it as follows:
I changed the constr routine to this:
subroutine constr( cell, stress, na, isa, amass, xa, fa, ntcon )
c *****************************************************************
use m_fdf_global, only: fdf_global_get
use parallel, only: IOnode
implicit none
integer na, isa(na), ntcon, ia1, ia2
double precision amass(na), cell(3,3), fa(3,na),
. stress, xa(3,na)
real af_magn
c Write here your problem-specific code.
call fdf_global_get(ia1,'MD.AFatom1',1)
call fdf_global_get(ia2,'MD.AFatom2',2)
call fdf_global_get(af_magn,'MD.AF',0.0)
if (IOnode) write(6,'(A)') ' '
if (IOnode) then
write(6,'(A, I4, I4)' )
. 'Bart: Artificial force added to atoms:', ia1, ia2
write(6, '(A)' ) 'Bart: Coordinates: '
write(6, '(A, I4, F12.5, F12.5, F12.5)' )
. 'Bart: ',ia1, xa(1, ia1), xa(2, ia1), xa(3,ia1)
write(6, '(A, I4, F12.5, F12.5, F12.5)' )
. 'Bart: ',ia2, xa(1, ia2), xa(2, ia2), xa(3,ia2)
write(6,'(A, F8.5)' ) 'Bart: Force magnitude: ', af_magn
write(6,'(A, F8.5)' ) 'Bart: Force magnitude: [nN] ',
. af_magn*41.1885
fa(3,ia1) = fa(3,ia1) + af_magn
fa(3,ia2) = fa(3,ia2) - af_magn
ntcon = 1
end
I needed to recompile SIESTA only once. Then I added MD.AF (force
magnitude), MD.AFatom1 and MD.AFatom2 (atoms to which the force is
applied) keywords to the fdf file. The lines with fdf_global_get in the
constr subroutine get this information from the fdf file. Then the
forces are modified - to the force acting on a particular atom resulting
from gradients and whatever else - the artificial force is added along
the z direction (defined by number 3 in fa(3,ia1), use 1 and 2 for x and
y respectively). You also get a few more lines about it in the output
file.
Feel free to use it and modify it according to your needs :)
Cheers,
Bartek
