module real_space_moire
     contains
     subroutine unit_cell_creation()
            implicit none
            double precision :: d1,d2,d
            double precision, dimension(2) :: a1, a2,  unit_cell_position
            integer :: unit_cell_ind1, unit_cell_ind2
            double precision :: cutoff_radius_inner, cutoff_radius_outer, theta
            double precision, dimension(2) :: tau_a, tau_b, pos1_long, pos2_long, U_j,  R_0, R_B1
            integer :: pre_lattice_size
            integer :: unit_cell_counter
            double precision, allocatable :: uc_pos(:,:)
            double precision, parameter :: pi=4.0d0*atan(1.0d0)

            cutoff_radius_inner=0.0d0
            cutoff_radius_outer=100.0d0

            d1=1.42

            R_0 = 0.0d0
            R_B1(1)=d1
            R_B1(2)=0.0d0

            !creation of lattice

            a1(1)=0.0d0
            a1(2)=d1*sqrt(3.0d0)

            a2(1)=(d1/2.0d0)*3.0d0
            a2(2)=(-d1/2.0d0)*sqrt(3.0d0)

            pre_lattice_size=1000+INT(cutoff_radius_outer)
            !cutoff_radius_inner = 0.00d0 !in Angst.
            !cutoff_radius_outer = 200.0d0 !in Angst.


            unit_cell_counter=0
            !lattice generation
            do unit_cell_ind1=-pre_lattice_size,pre_lattice_size
            do unit_cell_ind2=-pre_lattice_size,pre_lattice_size

                unit_cell_position(1) = unit_cell_ind1*a1(1)+unit_cell_ind2*a2(1)
                unit_cell_position(2) = unit_cell_ind1*a1(2)+unit_cell_ind2*a2(2)

                if ( sqrt(unit_cell_position(1)**2.0d0 + unit_cell_position(2)**2.0d0) .GE. cutoff_radius_inner .AND. &
                     sqrt(unit_cell_position(1)**2.0d0 + unit_cell_position(2)**2.0d0) .LE. cutoff_radius_outer   ) then
                unit_cell_counter=unit_cell_counter+1
                end if

            end do
            end do

            !write(*,*) 'unit_cell_counter', unit_cell_counter
            open (unit=10, file="atom_A.dat", recl=512)
            open (unit=11, file="atom_B.dat", recl=512)
            open (unit=12, file="atom_Ap.dat", recl=512)
            open (unit=13, file="atom_Bp.dat", recl=512)
            allocate( uc_pos(unit_cell_counter,2) )

            unit_cell_counter=1
            do unit_cell_ind1=-pre_lattice_size,pre_lattice_size
             do unit_cell_ind2=-pre_lattice_size,pre_lattice_size
                unit_cell_position(1) = unit_cell_ind1*a1(1)+unit_cell_ind2*a2(1)
                unit_cell_position(2) = unit_cell_ind1*a1(2)+unit_cell_ind2*a2(2)


                if ( sqrt(unit_cell_position(1)**2.0d0 + unit_cell_position(2)**2.0d0) .GE. cutoff_radius_inner .AND. &
                  sqrt(unit_cell_position(1)**2.0d0 + unit_cell_position(2)**2.0d0) .LE. cutoff_radius_outer   ) then

                    uc_pos(unit_cell_counter,1)=unit_cell_position(1)
                    uc_pos(unit_cell_counter,2)=unit_cell_position(2)

                    !atom A
                    write(10,*) uc_pos(unit_cell_counter,1)+R_0(1), uc_pos(unit_cell_counter,2)+R_0(2), 0.0d0
                    !atom B
                    write(11,*) uc_pos(unit_cell_counter,1)+R_B1(1), uc_pos(unit_cell_counter,2)+R_B1(2), 0.0d0

                    theta=0.3802512066 !60.0d0/360.0d0 * 2.0d0 *pi
                    !atom A'
                    write(12,*) cos(theta)*( uc_pos(unit_cell_counter,1)+R_0(1) -R_B1(1)) - sin(theta)*( uc_pos(unit_cell_counter,2)+R_0(2)- R_B1(2)), & ! x
                                sin(theta)*(uc_pos(unit_cell_counter,1) +R_0(1) -R_B1(1)) + cos(theta)*( uc_pos(unit_cell_counter,2)+R_0(2)- R_B1(2)), & ! y
                                -3.35d0 !z
                    !atom B'
                    write(13,*) cos(theta)*( uc_pos(unit_cell_counter,1) ) - sin(theta)* (uc_pos(unit_cell_counter,2)), & !x
                                sin(theta)*( uc_pos(unit_cell_counter,1) ) + cos(theta)* (uc_pos(unit_cell_counter,2)), & !y
                                -3.35d0 !z
                    unit_cell_counter=unit_cell_counter+1
                end if

              end do
            end do
            close(10)
            close(11)
            close(12)
            close(13)

            unit_cell_counter=unit_cell_counter-1
        end subroutine unit_cell_creation

    end module real_space_moire