Hi Roberto, Thank you for your help. I think you solved my issue. I missed, that the conjugate gradient algorithm I use is just an energy minimization and not a molecular dynamics simulation, thus doesn't follow the newton's laws of motion. So as far as I understood I should be able to safely ignore the mass in the constraint routine, as long as I use CG as minimization method.
Btw. I'm sorry that my last reply was a bit confusing about my main simulation. I don't want to apply that constraint to an isolated molecule, but to some nucleobases in a DNA. Thanks again for your help, Best wishes, Frank > Hi Frank, > > I've used the constraint routine many times and it ever worked fine for > me. > Your molecule is isolated, it should not rotate anyway becuase internal > forces possess no mechanical torque. I mean, such a constraint would be > superfluous. > I any case, if you still want to apply such a constraint, you should > find out how the atoms would move in such a virtual rotation, namely, > build up the rotation mode, then project the forces onto that mode, and > finally subtract the projection from the force vector. > Cheers, > > Roberto > > PS: Note that your're doing a static energy minimization, the masses play > no physical role here. > > On 08/13/2013 09:20 AM, Frank Maier wrote: >> 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 >>>> >>>> >>>> >> >
