Hi Roberto,
you're right, the setup is idiotic :-) but it's indeed what I want. It's
just a test setup to validate the constraint routine and it nicely shows
that it doesn't work.
For my proper simulation I want to block the rotation of the molecule
along one axis with the constraint routine but still allow a
translation. For this I calculated the sum of all forces acting on the
molecule along the z-axis and distributed this force over all atoms,
identical to how it gets done in the example in the SIESTA manual. But
then I noticed, that the hydrogen atoms in the molecule barely move at
all while the heavier atoms move away.
Thus I constructed this artificial system to find the failure. But now
I'm stuck.
Including the mass makes sense to me, but it doesn't work.
Neglecting the mass works, but doesn't make sense to me.
Best wishes,
Frank
On 08/13/2013 02:11 PM, Roberto Pasianot wrote:
Hi Frank,
What do you want to do ?. A net force applied to an isolated
body causes the latter to never stop moving, it doesn't matter
how the force is distributed. This is surely not what you
mean.
Regards,
Roberto
On 08/13/2013 05: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
