#! /usr/bin/env python

import sys
import string
import math
from random import gauss
from numpy import linalg, dot
import numpy as np

POSCAR = open('POSCAR', 'r').readlines()

N=38
g=open('cfg_lammps.RuN.config','w')
g.write('Written by Genlammps_AIMD.py, original comment line\n\n')
g.write('   {0:.0f}  atoms\n'.format(N))
g.write('   3  atom types\n\n')
g.write( "  0.0000000000000000 8.6930937287162351 xlo xhi\n" )
g.write( "  0.0000000000000000 7.5284400065474486 ylo yhi\n" )
g.write( "  0.0000000000000000 22.2281052000000017 zlo zhi\n" )
g.write( "  4.3465468643581175 0.0000000000000000 0.0000000000000000 xy xz yz\n\n" )
g.write('Atoms\n\n')
POS = [[ 0.0, 0.0, 7.5], [1.279978, 0.0, 7.5]]
VEL = [[ 0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
for k in range(0,N):
	if k < 36:
		g.write( " {:4d}  3  {:16.12f}  {:16.12f}  {:16.12f}\n".format( k+1, float(POSCAR[k+9].split()[0]), float(POSCAR[k+9].split()[1]), float(POSCAR[k+9].split()[2]) ) )
	elif k == 36:
		g.write( " {:4d}  1  {:16.12f}  {:16.12f}  {:16.12f}\n".format( k+1, float(POS[0][0]), float(POS[0][1]), float(POS[0][2]) ) )
	else:
		g.write( " {:4d}  2  {:16.12f}  {:16.12f}  {:16.12f}\n".format( k+1, float(POS[1][0]), float(POS[1][1]), float(POS[1][2]) ) )

g.write('\n')
g.write('Velocities\n\n')
for k in range(0,N-2):
	g.write( " {:4d}  {:16.12f}  {:16.12f}  {:16.12f}\n".format( k+1, float(POSCAR[k+46].split()[0])*1000, float(POSCAR[k+46].split()[1])*1000, float(POSCAR[k+46].split()[2])*1000 ) )
g.write( " {:4d}  {:16.12f}  {:16.12f}  {:16.12f}\n".format( N-1, float(VEL[0][0]), float(VEL[0][1]), float(VEL[0][2]) ) )
g.write( " {:4d}  {:16.12f}  {:16.12f}  {:16.12f}\n".format( N, float(VEL[1][0]), float(VEL[1][1]), float(VEL[1][2]) ) )
g.close()