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
