Thank you for your answer Valentin,
but sadly I don't fully understand it. In the constraint routine I
overwrite the forces acting on the atoms along the z-axis with my custom
defined ones. What other forces are there, or what do you mean with
external forces?
I also don't understand, even in the presence of other, to me unknown,
forces, why it behaves correct (but in my opinion isn't physically
correct any longer) if I apply the identical force to each atom
regardless of their mass and it doesn't behave correct if I apply a
weighted (depending on their mass) force (which is physically correct in
my opinion) to each atom. In the case of the weighted force, all atoms
with the same mass move homogeneously.
Best wishes,
Frank
On 08/13/2013 09:54 AM, Valentin Karasiev wrote:
Frank,
do not forget, the external force you apply is not the only force
acting on atoms. In the case of bigger external force, the interatomic
forces may be neglected. In the case of smaller external force,
interatomic forces should be taken into account. It seems you are
observing namely that in your calculations.
Valentin.
On 08/13/2013 04:49 AM, Frank Maier wrote:
Hello,
I try to add a constraint to my system using the constraint routine
in Fortran using the latest SIESTA Trunk 433.
Sadly I encountered an issue which I was able to break down to a
simple test system:
I want to apply a force of 1eV/Ang along the z-axis homogeneously to
the DNA base Thymine which consists of 15 atoms in total. In order to
do this, I first calculate the sum of the mass of all atoms and then
apply a weighted force along the z-axis to each atom depending on
their mass relative to the whole mass of the molecule (the code is
attached at the end of the mail).
If I do this, the hydrogen atoms barely move, oxygen atoms move the
most, followed by nitrogen and carbon. So depending on their mass,
they move differently.
If I apply 1/15 eV/Ang to each atom however, they move equivalently,
which shouldn't happen in my opinion. If I apply a force of 1N to a
mass of 1kg it should have a different impact compared to a force of
1N to a mass of 16kg.
In principle, for my real simulation, i want to preserve the relative
height of atoms in a molecule, identical to how it gets done in the
Siesta manual example:
fz = fa(3,1) + fa(3,2)
fa(3,1) = fz * amass(1)/(amass(1)+amass(2))
fa(3,2) = fz * amass(2)/(amass(1)+amass(2))
But this suffers the same problem, the atoms just move differently,
depending on their mass. So exactly the opposite to how it should
happen.
Can someone help me.
Below is the constraint routine I use attached and also the
simulation file so you can understand what I've done.
Thank you for your help,
Frank Maier
########### CONSTRAINT ROUTINE ###########
implicit none
integer na, isa(na), ntcon
double precision amass(na), cell(3,3), fa(3,na),
. stress(3,3), xa(3,na)
double precision force_z
double precision mass_sum
integer counter
c Write here your problem-specific code
! sum all masses
mass_sum = 0
counter = 1
do while (counter <= 15)
mass_sum = mass_sum + amass(counter)
counter = counter + 1
end do
! drag it along the z axis by setting an average force to all atoms
force_z = 0.038895
counter = 1
do while (counter <= 15)
fa(3, counter) = force_z * amass(counter)/mass_sum
counter = counter + 1
end do
! number of constraints
ntcon = 1
end
########### SIESTA INPUT FILE ###########
SystemName simulation
SystemLabel simulation
NumberOfAtoms 15
NumberOfSpecies 4
%block Chemical_Species_label
1 6 C
2 1 H
3 7 N
4 8 O
%endblock Chemical_Species_label
LatticeConstant 1 Ang
%block LatticeVectors
20.0 0.0 0.0
0.0 20.0 0.0
0.0 0.0 20.0
%endblock LatticeVectors
AtomicCoordinatesFormat NotScaledCartesianAng
%block AtomicCoordinatesAndAtomicSpecies
-0.055415 4.678721 0.046505 4 1 O
0.056186 0.053516 0.044754 4 2 O
-2.040067 3.505861 -0.185110 3 3 N
-0.021091 2.357257 0.016693 3 4 N
-0.646655 3.608724 -0.032440 1 5 C
-0.631420 1.071320 -0.034209 1 6 C
-2.096698 1.092563 -0.183565 1 7 C
-2.817098 -0.233357 -0.250395 1 8 C
-2.727280 2.298413 -0.253138 1 9 C
0.990177 2.373717 0.118928 2 10 H
-2.534782 4.389011 -0.226380 2 11 H
-3.804217 2.386022 -0.365223 2 12 H
-3.898566 -0.088334 -0.347702 2 13 H
-2.456966 -0.824419 -1.102264 2 14 H
-2.614350 -0.824912 0.651503 2 15 H
%endblock AtomicCoordinatesAndAtomicSpecies
%block GeometryConstraints
routine constr
%endblock GeometryConstraints
PAO.BasisSize TZP
PAO.BasisType split
PAO.EnergyShift 10.0 meV
PAO.SplitNorm 0.20
PAO.SplitNormH 0.50
PAO.SoftDefault true
PAO.SoftInnerRadius 0.9
PAO. SoftPotential 40.0 Ry
XC.functional GGA
XC.authors PBE
SpinPolarized false
SolutionMethod diagon
ElectronicTemperature 300 K
NeglNonOverlapInt false
DM.MixingWeight 0.1
DM.Tolerance 1.d-4
DM.NumberPulay 4
MeshCutoff 250. Ry
MD.TypeOfRun CG
MD.NumCGsteps 500
MD.MaxCGDispl 0.2 Bohr
MD.MaxForceTol 0.02 eV/Ang
