!/***************************************************************************
! *   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 fcc lattice for a squared simulation box
!
!*****************************************************

program makefcc
implicit none

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

write(*,*) ' * * * Generate a FCC 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')


!****** Lattice Vectors ************
      v1(1)=0.d0
      v1(2)=0.5d0
      v1(3)=0.5d0
     
      v2(1)=0.5d0
      v2(2)=0.d0
      v2(3)=0.5d0
      
      v3(1)=0.5d0
      v3(2)=0.5d0
      v3(3)=0.d0
!****** Lattice Vectors ************

! This condition is valid only for squared simulation box

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

      walk=1

      delta=LBox*(4.d0/dble(nion))**(1.d0/3.d0)
     
      L2=LBox/2.d0

      do i = -nion/2,nion/2
       do j = -nion/2,nion/2
        do k = -nion/2,nion/2
         rt(1)=delta*(v1(1)*i+v2(1)*j+v3(1)*k) 
         rt(2)=delta*(v1(2)*i+v2(2)*j+v3(2)*k) 
         rt(3)=delta*(v1(3)*i+v2(3)*j+v3(3)*k) 
         if(rt(1).lt.L2.and.rt(2).lt.L2.and.rt(3).lt.L2&
    .and.rt(1).ge.-L2.and.rt(2).ge.-L2.and.rt(3).ge.-L2) then
           rion(1,walk)=rt(1)
           rion(2,walk)=rt(2)
           rion(3,walk)=rt(3)
           walk=walk+1
        endif
        enddo
       enddo
      enddo
      
write(11,*) nion
write(11,*) '# Fcc Lattice ' 
do i=1,nion
        write(11,*) ' C ',(rion(j,i),j=1,3)
enddo

      end