!/***************************************************************************
! *   Copyright (C) 2005 by Claudio Attaccalite                             *
! *   claudio@freescience.info                                              *
! *                                                                         *
! *   This program is free software; you can redistribute it and/or modify  *
! *   it under the terms of the GNU General Public License as published by  *
! *   the Free Software Foundation; either version 2 of the License, or     *
! *   (at your option) any later version.                                   *
! *                                                                         *
! *   This program is distributed in the hope that it will be useful,       *
! *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
! *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
! *   GNU General Public License for more details.                          *
! *                                                                         *
! *   You should have received a copy of the GNU General Public License     *
! *   along with this program; if not, write to the                         *
! *   Free Software Foundation, Inc.,                                       *
! *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
! ***************************************************************************/

!*****************************************************
!
! Generate an bcc lattice for a squared simulation box
!
!*****************************************************

program makebcc
implicit none

integer nion,N,walk,i,j,k
integer done,ct
real(8), DIMENSION(:,:), ALLOCATABLE :: rion
real(8) maxdist,dist,rs,LBox,delta,L2

write(*,*) ' * * * Generate a BCC Lattice for a squared simulation cell * * * '
write(*,*)
write(*,'(" Rs = ")',ADVANCE='NO')
read(*,*) rs
write(*,'(" Number of ions: ")',ADVANCE='NO')
read(*,*) nion

LBox = (3.14159265358979323846*nion*4.d0/3.d0)**(1.d0/3.d0)*rs
write(*,*) ' LBox = ',LBox

ALLOCATE(rion(3,nion))
     
!********* Output file ***********
open(unit=11,file="output.xyz",form='formatted',status='unknown')


! This condition is valid only for squared simulation box

  N=int((nion/2.d0)**(1.d0/3.d0)+0.5d0)
  if(N*N*N.ne.nion/2.d0) then
    write(*,*) " Error! Not valid number of ions for a squared simulation box "
    write(*,*) " Number of ions = ",nion," choose=  ",N*N*N*2 
    stop
  endif

      walk=1
      delta=LBox/dble(N)
  
      do k = 0,N-1
       do j = 0,N-1
        do i = 0,N-1
         rion(1,2*walk)= -0.5*LBox+delta*i+delta/4.d0;
         rion(2,2*walk)= -0.5*LBox+delta*j+delta/4.d0;
         rion(3,2*walk)= -0.5*LBox+delta*k+delta/4.d0;
         rion(1,2*walk-1)= -0.5*LBox+delta*i+delta*3.0/4.0;
         rion(2,2*walk-1)= -0.5*LBox+delta*j+delta*3.0/4.0;
         rion(3,2*walk-1)= -0.5*LBox+delta*k+delta*3.0/4.0;
         walk=walk+1;
          enddo
        enddo
      enddo
      
write(11,*) nion
write(11,*) '# Bcc Lattice ' 
do i=1,nion
        write(11,*) ' C ',(rion(j,i),j=1,3)
enddo

      end