import numpy as np
from ase import Atoms
from ase.io import read
from ase.lattice.spacegroup import crystal
from ase.optimize import BFGS, BFGSLineSearch
from ase.constraints import FixAtoms, UnitCellFilter
from gpaw import GPAW, PW, MixerSum, MixerDif, FermiDirac
from gpaw.occupations import FermiDirac

#From gs.gpw
pos = [[  3.14635755e+00,   1.81655038e+00,   8.87569021e+00],
       [  1.08133308e-15,   3.63310076e+00,   8.87566200e+00],
       [  3.97192780e+00,  -2.06299436e-04,   7.49522829e+00],
       [  1.16021499e+00,   2.00996392e+00,   7.49522829e+00],
       [ -1.98578524e+00,   3.43989352e+00,   7.49522829e+00],
       [  2.32078745e+00,  -1.95870631e-04,   1.02561218e+01],
       [  1.98579420e+00,   3.43988819e+00,   1.02561218e+01],
       [ -1.16022409e+00,   2.00995882e+00,   1.02561218e+01]]

cellp = [[  6.29271511e+00,  -1.57210150e-19,   0.00000000e+00],
      [ -3.14635755e+00,   5.44965114e+00,   0.00000000e+00],
       [ -7.91811112e-19,   0.00000000e+00,   1.77513506e+01]]

a = Atoms('Fe2Cl6', cell=cellp, positions = pos, pbc = [True, True, False] )

magmom = [ 1,  1,  0,  0,  0,  0,  0,  0]


a.set_initial_magnetic_moments(magmom)

for i in range(200,901, 50):
	calc = GPAW(mode=PW(i),
		    xc='PBE',
		    #mixer=MixerDif(),
		    occupations=FermiDirac(0.1),
		    parallel={'domain': 1, 'band': 1},
		    kpts=(8,8,1),
		    txt='conv_pw_FeCl3%d.txt' % i)
	a.set_calculator(calc)
	a.get_potential_energy()

	uf = UnitCellFilter(a)
	#opt = BFGSLineSearch(uf)
	opt = BFGS(uf)
	opt.run(fmax=0.01)
	calc.write('conv_pw_FeCl3%d.gpw' % i)

#for k in range(2, 13, 2):
#	calc = GPAW(mode=PW(600),
#		    xc='PBE',
#		    #mixer=MixerDif(),
#		    occupations=FermiDirac(0.1),
#		    parallel={'domain': 1, 'band': 1},
#		    kpts=(k,k,1),
#		    txt='conv_k_FeCl3%d.txt' % k)
#	a.set_calculator(calc)
#	a.get_potential_energy()
#
#	uf = UnitCellFilter(a)
	#opt = BFGSLineSearch(uf)
#	opt = BFGS(uf)
#	opt.run(fmax=0.01)
#	calc.write('conv_k_FeCl3%d.gpw' % k)