Michael,
What are the forces calculated by SIESTA during the run? You can find
them in the output file. If you are interested in implementation of the
external forces in SIESTA, you can try my code (Szyja, B.M., Pidko,
E.A., Groote, R., Hensen, E.J.M. & Sijbesma, R.P. (2011). DFT study on
mechanochemical bond breaking in COGEF and molecular dynamics
simulations. Procedia Computer Science, 4, 1167-1176.) It has this
advantage, that it looks for the parameters in the fdf file, and you
don't need to recompile SIESTA every time you want to change the
magnitude of force. It is very simple but worked well for me, although
my goal was not to fix the position of atoms.
Cheers,
Bartek
On Sat, 2012-05-12 at 21:50 +0200, Michael Shin wrote:
> Rurali,
> I llokedinto your discussion and it helped us. Then I tried Silican as
> an example.
> I have two atoms and I wanted not to relax both the atoms and modified
> the constr.f file andthen recomplie it, but it didnt wirk.
> Then I thought lets add constrain on two atoms and I added the
> following line in the constr.f file.
> fa(3,1)=0.0d0
> fa(3,2)=0.0d0
>
> Which means that folrce along z axis on atom 1 and two will be zero.
> When I recompiled the code and run it, thn it relaxed both the atoms
> in the unit cell.
> I started from the following lines
>
> %block AtomicCoordinatesAndAtomicSpecies
> 0. 0. 0.3 1 Si 1
> 0.25 0.25 0.26 1 Si 2
> %endblock AtomicCoordinatesAndAtomicSpecies
>
> And got the Following relaxed cocordines.
>
> outcoor: Relaxed atomic coordinates (fractional):
> 0.00054664 0.00054723 0.15300834 1 1 Si
> 0.25058586 0.25058642 0.40284097 1 2 Si
>
> In principles,if I am right, it should not relax these two atoms, but
> it does.
> Can you help me how to handle it. Actually my system has about 100
> atoms and I dont need to relax all the atoms, just a few atoms near
> the surface. So I am learning this method and tries Silicon.
>
> Mic
>
>
> --- On Fri, 1/21/11, Riccardo Rurali <[email protected]> wrote:
>
>
> From: Riccardo Rurali <[email protected]>
> Subject: Re: [SIESTA-L] For help: How to pply a constant
> initial velocity or external force on some atoms of a system.
> To: [email protected]
> Date: Friday, January 21, 2011, 6:36 AM
>
> Hi there.
>
> On 1/20/11 5:29 PM, Li, Shunfang wrote:
> > My questions are the following:
> > 1)Can we apply a constant initial velocity on some atoms of
> a system, in
> > a MD simulations performed by Siesta?
> > I know, the atoms are assigned random velocities drawn from
> the
> > Maxwell-Bolzmann distribution with the corresponding
> temperature in a MD
> > simulation, however, you know, the constraint of zero center
> of mass
> > velocity is imposed.
>
> If you want to apply a constant *initial* velocity this is
> easily done with a little trick.
> You run a fake calculation at your interested target
> temperature. I call it "fake", because you need this
> calculation only to have Siesta write the XV file. You don't
> even need to converge the electronic density.
>
> The XV file will contain the atomic positions (your initial
> coordinates if you don't relax) and the velocities sampled
> form the Maxwell-Boltzmann distribution, as you said.
>
> You can edit that file and define the velocity you like for
> the atoms(s) you like.
> You are now ready to run your calculation, making sure you
> specify the option MD.UseSaveXV true. In this way the initial
> positions and velocities are read from the XV file.
>
> This should do what you want.
>
> > 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
>
> Good luck,
> Riccardo
>
>
> --
> Riccardo Rurali
> Institut de Ciència de Materials de Barcelona (ICMAB)
> Consejo Superior de Investigaciones CientÃficas (CSIC)
> Campus de Bellaterra
> 08193 Bellaterra (Barcelona)
> Spain
>
> tel.: +34 93 5801853 ext. 347
> e-mail: [email protected]
> http://www.icmab.es/dmmis/leem/
>
> Man, the dope's that there's still hope
>
>
!
! This file is part of the SIESTA package.
!
! Copyright (c) Fundacion General Universidad Autonoma de Madrid:
! E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal
! and J.M.Soler, 1996- .
!
! Use of this software constitutes agreement with the full conditions
! given in the SIESTA license, as signed by all legitimate users.
!
c $Id: constr.f,v 1.6 2003/06/23 09:46:16 ordejon Exp $
subroutine constr( cell, stress, na, isa, amass, xa, fa, ntcon )
c *****************************************************************
c User-written routine to implement specific geometric constraints,
c by orthogonalizing the forces and stress to undesired changes.
c Arguments:
c real*8 cell(3,3) : input lattice vectors (Bohr)
c integer na : input number of atoms
c integer isa(na) : input species indexes
c real*8 amass(na) : input atomic masses
c real*8 xa(3,na) : input atomic cartesian coordinates (Bohr)
c real*8 stress( 3,3) : input/output stress tensor (Ry/Bohr**3)
c real*8 fa(3,na) : input/output atomic forces (Ry/Bohr)
c integer ntcon : total number of positions constr. imposed
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
write(6,'(A)') 'Bart: Original forces: [nN] '
write(6,'(A,I4,F12.5,F12.5,F12.5)') 'Bart: ',
. ia1, fa(1,ia1)*41.1885, fa(2,ia1)*41.1885, fa(3,ia1)*41.1885
write(6,'(A,I4,F12.5,F12.5,F12.5)') 'Bart: ',
. ia2, fa(1,ia2)*41.1885, fa(2,ia2)*41.1885, fa(3,ia2)*41.1885
endif
fa(3,ia1) = fa(3,ia1) + af_magn
fa(3,ia2) = fa(3,ia2) - af_magn
if (IOnode) then
write(6,'(A)') 'Bart: Updated forces: [nN] '
write(6,'(A,I4,F12.5,F12.5,F12.5)') 'Bart: ',
. ia1, fa(1,ia1)*41.1885, fa(2,ia1)*41.1885, fa(3,ia1)*41.1885
write(6,'(A,I4,F12.5,F12.5,F12.5)') 'Bart: ',
. ia2, fa(1,ia2)*41.1885, fa(2,ia2)*41.1885, fa(3,ia2)*41.1885
endif
ntcon = 1
end
