Hi all
I am trying to optimize a crystal with some restrictions, such alpha,
beta, and gamma angles are constants. Moreover, I also want to fix the
unit-cell side lengths to their initial values for b and c. This is my
file constr.f:
subroutine constr( cell, na, isa, amass, xa, stress, fa, ntcon )
c *****************************************************************
c User-written routine to implement specific geometric constraints,
c by orthogonalizing the forces and stress to undesired changes.
c Arguments:
c real*8 cell(3,3) : input lattice vectors (Bohr)
c integer na : input number of atoms
c integer isa(na) : input species indexes
c real*8 amass(na) : input atomic masses
c real*8 xa(3,na) : input atomic cartesian coordinates (Bohr)
c real*8 stress( 3,3) : input/output stress tensor (Ry/Bohr**3)
c real*8 fa(3,na) : input/output atomic forces (Ry/Bohr)
c integer ntcon : total number of positions constr. imposed
c *****************************************************************
implicit none
integer na, isa(na), ntcon, ia
double precision amass(na), cell(3,3), fa(3,na),
. stress(3,3), xa(3,na)
double precision fav, stressav
c Write here your problem-specific code.
stress(1,2) = 0.0d0
stress(1,3) = 0.0d0
stress(2,1) = 0.0d0
stress(2,2) = 0.0d0
stress(2,3) = 0.0d0
stress(3,1) = 0.0d0
stress(3,2) = 0.0d0
stress(3,3) = 0.0d0
end
I compile the programm and finally, I add the block label:
%block GeometryConstraints
routine constr
%endblock
I show the end part of output file:
siesta: ==============================
Begin CG move = 0
==============================
superc: Internal auxiliary supercell: 2 x 2 x 1 = 4
superc: Number of atoms, orbitals, and projectors: 400 3408 4832
outcell: Unit cell vectors (Ang):
7.872628 0.000000 0.000000
0.000000 7.800000 0.000000
0.000000 0.000000 17.000000
outcell: Cell vector modules (Ang) : 7.872628 7.800000 17.000000
outcell: Cell angles (23,13,12) (deg): 90.0000 90.0000 90.0000
outcell: Cell volume (Ang**3) : 1043.9105
InitMesh: MESH = 96 x 96 x 216 = 1990656
InitMesh: Mesh cutoff (required, used) = 400.000 410.965 Ry
* Maximum dynamic memory allocated = 630 MB
stepf: Fermi-Dirac step function
siesta: Program's energy decomposition (eV):
siesta: Eions = 13406.371149
siesta: Ena = 3017.931671
siesta: Ekin = 4936.047104
siesta: Enl = -22.576055
siesta: DEna = -0.000044
siesta: DUscf = 0.000000
siesta: DUext = 0.000000
siesta: Exc = -2602.816171
siesta: eta*DQ = 0.000000
siesta: Emadel = 0.000000
siesta: Ekinion = 0.000000
siesta: Eharris = -8216.934182
siesta: Etot = -8077.784645
siesta: FreeEng = -8077.784645
siesta: iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV)
siesta: 1 -8216.9342 -8077.7846 -8077.7846 1.9478 -4.8575
timer: Routine,Calls,Time,% = IterSCF 1 225.838 95.38
elaps: Routine,Calls,Wall,% = IterSCF 1 225.937 95.38
siesta: 2 -8225.2616 -8186.5448 -8186.5448 0.2923 -2.8212
siesta: 3 -8203.2830 -8183.6343 -8183.6343 0.1624 -3.3919
siesta: 4 -8201.6096 -8194.8313 -8194.8313 0.0766 -3.4613
siesta: 5 -8201.5367 -8198.0680 -8198.0680 0.0457 -3.4062
siesta: 6 -8201.5311 -8200.7177 -8200.7177 0.0133 -3.4112
siesta: 7 -8201.5274 -8201.2270 -8201.2270 0.0029 -3.4388
siesta: 8 -8201.5271 -8201.3701 -8201.3701 0.0011 -3.4355
siesta: 9 -8201.5271 -8201.4526 -8201.4526 0.0005 -3.4331
siesta: 10 -8201.5271 -8201.5118 -8201.5118 0.0002 -3.4336
siesta: 11 -8201.5270 -8201.5294 -8201.5294 0.0001 -3.4338
siesta: E_KS(eV) = -8201.5299
siesta: E_KS - E_eggbox = -8201.5299
Which is the problem?
Thanks in advance
Gregorio Garcia
