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
