updated

src/George_Marsaglia.f90

@@ -0,0 +1,39 @@
+! Ether, a Monte Carlo simulation program, impowers the users to study 
+! the thermodynamical properties of spins arranged in any complex 
+! lattice network.
+
+! Copyright (C) 2021  Mukesh Kumar Sharma ([email protected])
+
+! 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, see https://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html
+
+	subroutine George_Marsaglia(spinx, spiny, spinz)
+
+		implicit none
+
+		real(8), intent(out) :: spinx, spiny, spinz
+		real(8) :: v1, v2, S_val, rn
+
+7		call get_random_num(-1d0, +1d0, rn)
+		v1 = rn
+		call get_random_num(-1d0, +1d0, rn)
+		v2 = rn
+		S_val = v1**2 + v2**2
+
+		if (S_val.ge.1) goto 7
+
+		spinx = 2*v1*sqrt(1-S_val)
+		spiny = 2*v2*sqrt(1-S_val)
+		spinz = 1-2*S_val
+
+	end subroutine George_Marsaglia

src/Hamiltonian.f90

@@ -0,0 +1,138 @@
+! Ether, a Monte Carlo simulation program, impowers the users to study 
+! the thermodynamical properties of spins arranged in any complex 
+! lattice network.
+
+! Copyright (C) 2021  Mukesh Kumar Sharma ([email protected])
+
+! 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, see https://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html
+
+	subroutine Hamiltonian(ih, jh, kh, lh, Si_center, total_ene)
+
+		use init
+
+		implicit none
+
+		integer, intent(in) :: ih, jh, kh, lh
+		integer :: ionID, nbdi, inbr, cell
+
+		real(8) :: sia_value(3)
+		real(8), intent(in) :: Si_center(3)	! central ion
+		real(8), intent(out) :: total_ene
+
+		total_ene = 0
+
+		! Central ion's species ID
+		ionID = int(ion(4, ih, jh, kh, lh))
+
+		call JSiSj(ih, jh, kh, lh, Si_center, ionID, total_ene)
+
+		if(Zeeman) call gmbSH(Si_center, g_factor, mb, H, total_ene)
+
+		if(single_ion_anisotropy.and..not.Ising) then
+			sia_value(1:3) = sia_vec(1:3, ionID)
+			call siaS2(Si_center, sia_value, total_ene)
+		end if
+
+	contains
+
+	! JSiSj	
+	subroutine JSiSj(i, j, k, l, Si, central_ion_ID, total_energy)
+
+		implicit none
+
+		integer, intent(in) :: i, j, k, l, central_ion_ID
+		integer :: total_nbr, posx, posy, posz, ID_num, ion_ID
+
+		real(8), intent(in) :: Si(3)
+		real(8) :: Sj(3), SiSj(3), Jij(3), eout
+		real(8), intent(inout) :: total_energy
+
+		ion_ID = int(ion(0, i, j, k, l))
+		! Over distinct bonds
+		do nbdi = 1, no_of_nbd
+
+			! For nbdi bond total connecting neighbours to the central ion [ID = ion(0, i, j, k, l)] is
+			total_nbr = nn(nbdi, ion_ID, 0, 0)
+
+			! no. of similar nbd for ith distinct bond (nbdi)
+			do inbr = 1, total_nbr
+
+				! nbr's cell position
+				posx = nn(nbdi, ion_ID, inbr, 1)
+				posy = nn(nbdi, ion_ID, inbr, 2)
+				posz = nn(nbdi, ion_ID, inbr, 3)
+
+				! nbr's ID in the cell
+				cell = nn(nbdi, ion_ID, inbr, 4)
+				ID_num = int(ion(4, posx, posy, posz, cell))
+				if(ID_num.eq.0) then
+					go to 4
+				end if
+
+				Sj(1:3) = ion(1:3, posx, posy, posz, cell)
+
+				!Si.Sj
+				SiSj = Si*Sj	!point-wise multiplication
+
+				!Jij term
+				Jij = j_exc(nbdi, central_ion_ID, ID_num, 1:3)
+
+				!Jij.Si.Sj term
+				eout = dot_product(SiSj, Jij)
+
+				total_energy = total_energy + eout
+
+4				continue
+			end do
+		end do
+
+	end subroutine JSiSj
+
+	!gmbSH	
+	subroutine gmbSH(Si, g, mu, H, total_energy)
+
+		implicit none
+
+		real(8) :: eout
+		real(8), intent(in) :: Si(3), g, mu, H(3)
+		real(8), intent(inout) :: total_energy
+
+		! energy due to magnetic field
+		eout = -(g*mu*dot_product(Si, H))
+
+		total_energy = total_energy + eout
+
+	end subroutine gmbSH
+
+	!SINGLE ION ANISOTRPY (SIA)
+	subroutine siaS2(Si, sia_val, total_energy)
+
+		implicit none
+
+		real(8) :: Si2(3), eout
+		real(8), intent(in) :: Si(3), sia_val(3)
+		real(8), intent(inout) :: total_energy
+
+		! Si**2
+		Si2 = Si**2
+
+		! energy due to single ion anisiatropy
+		eout = dot_product(Si2, sia_val)			
+
+		total_energy = total_energy + eout
+
+	end subroutine siaS2
+
+	end subroutine Hamiltonian
+

src/Makefile

@@ -0,0 +1,3 @@
+ROOTDIR=..
+
+include ../Makefile

src/Makefile.ether

@@ -0,0 +1,125 @@
+# Avoid making changes below this line
+
+include ../make.sys
+
+OBJS  = init.o get_random_num.o read_input.o read_structure.o get_tot_species.o boundary_condition.o lu.o j_values.o \
+	parameters.o George_Marsaglia.o update_spin_details.o get_sia_values.o generate_supercell.o write_initial_conf.o \
+	getting_nbd.o write_nbd.o startup.o generate_bubble_indices.o Hamiltonian.o get_tot_energy.o \
+	get_tot_magnetisation.o fresh_spins.o allocate_observables.o evaluate_observables.o \
+	Monte_Carlo.o mc.o write_gss.o get_moment_vectors.o write_output_files.o write_spins_at_K.o zeroes.o \
+	get_ovrr_vec.o overrelaxation.o graph.o
+
+ether: ../ether
+
+../ether: $(OBJS) ../src/ether.f90
+	$(f90comp) $(switch) $(OBJS) ../src/ether.f90 -o ../ether
+
+clean:
+	rm -f *.o *.mod *.MOD *.obj
+
+init.o: ../src/init.f90
+	$(f90comp) -c $(switch) ../src/init.f90
+
+startup.o: ../src/startup.f90
+	$(f90comp) -c $(switch) ../src/startup.f90
+
+get_random_num.o: ../src/get_random_num.f90 init.o
+	$(f90comp) -c $(switch) ../src/get_random_num.f90
+
+read_input.o: ../src/read_input.f90 init.o
+	$(f90comp) -c $(switch) ../src/read_input.f90
+
+read_structure.o: ../src/read_structure.f90
+	$(f90comp) -c $(switch) ../src/read_structure.f90
+
+get_tot_species.o: ../src/get_tot_species.f90
+	$(f90comp) -c $(switch) ../src/get_tot_species.f90
+
+boundary_condition.o: ../src/boundary_condition.f90
+	$(f90comp) -c $(switch) ../src/boundary_condition.f90
+
+lu.o: ../src/lu.f90
+	$(f90comp) -c $(switch) ../src/lu.f90
+
+j_values.o: ../src/j_values.f90 init.o
+	$(f90comp) -c $(switch) ../src/j_values.f90
+
+parameters.o: ../src/parameters.f90 init.o
+	$(f90comp) -c $(switch) ../src/parameters.f90
+
+George_Marsaglia.o: ../src/George_Marsaglia.f90
+	$(f90comp) -c $(switch) ../src/George_Marsaglia.f90
+
+update_spin_details.o: ../src/update_spin_details.f90
+	$(f90comp) -c $(switch) ../src/update_spin_details.f90
+
+get_sia_values.o: ../src/get_sia_values.f90 init.o
+	$(f90comp) -c $(switch) ../src/get_sia_values.f90
+
+generate_supercell.o: ../src/generate_supercell.f90 init.o
+	$(f90comp) -c $(switch) ../src/generate_supercell.f90
+
+write_initial_conf.o: ../src/write_initial_conf.f90 init.o
+	$(f90comp) -c $(switch) ../src/write_initial_conf.f90
+
+getting_nbd.o: ../src/getting_nbd.f90 init.o
+	$(f90comp) -c $(switch) ../src/getting_nbd.f90
+
+write_nbd.o: ../src/write_nbd.f90 init.o
+	$(f90comp) -c $(switch) ../src/write_nbd.f90
+
+generate_bubble_indices.o: ../src/generate_bubble_indices.f90 init.o
+	$(f90comp) -c $(switch) ../src/generate_bubble_indices.f90
+
+Hamiltonian.o: ../src/Hamiltonian.f90 init.o
+	$(f90comp) -c $(switch) ../src/Hamiltonian.f90
+
+get_tot_energy.o: ../src/get_tot_energy.f90 init.o
+	$(f90comp) -c $(switch) ../src/get_tot_energy.f90
+
+get_tot_magnetisation.o: ../src/get_tot_magnetisation.f90 init.o
+	$(f90comp) -c $(switch) ../src/get_tot_magnetisation.f90
+
+fresh_spins.o: ../src/fresh_spins.f90 init.o
+	$(f90comp) -c $(switch) ../src/fresh_spins.f90
+
+allocate_observables.o: ../src/allocate_observables.f90 init.o
+	$(f90comp) -c $(switch) ../src/allocate_observables.f90
+
+evaluate_observables.o: ../src/evaluate_observables.f90 init.o
+	$(f90comp) -c $(switch) ../src/evaluate_observables.f90
+
+Monte_Carlo.o: ../src/Monte_Carlo.f90 init.o
+	$(f90comp) -c $(switch) ../src/Monte_Carlo.f90
+
+mc.o: ../src/mc.f90 init.o
+	$(f90comp) -c $(switch) ../src/mc.f90
+
+write_gss.o: ../src/write_gss.f90 init.o
+	$(f90comp) -c $(switch) ../src/write_gss.f90
+
+write_inital_conf.o: ../src/write_inital_conf.f90 init.o
+	$(f90comp) -c $(switch) ../src/write_inital_conf.f90
+
+get_moment_vectors.o: ../src/get_moment_vectors.f90 init.o
+	$(f90comp) -c $(switch) ../src/get_moment_vectors.f90
+
+write_output_files.o: ../src/write_output_files.f90 init.o
+	$(f90comp) -c $(switch) ../src/write_output_files.f90
+
+write_spins_at_K.o: ../src/write_spins_at_K.f90 init.o
+	$(f90comp) -c $(switch) ../src/write_spins_at_K.f90
+
+zeroes.o: ../src/zeroes.f90 init.o
+	$(f90comp) -c $(switch) ../src/zeroes.f90
+
+get_ovrr_vec.o: ../src/get_ovrr_vec.f90 init.o
+	$(f90comp) -c $(switch) ../src/get_ovrr_vec.f90
+
+overrelaxation.o: ../src/overrelaxation.f90 init.o
+	$(f90comp) -c $(switch) ../src/overrelaxation.f90
+
+graph.o: ../src/graph.f90 init.o
+	$(f90comp) -c $(switch) ../src/graph.f90
+
+.PHONY: ether clean

src/Monte_Carlo.f90

@@ -0,0 +1,85 @@
+! Ether, a Monte Carlo simulation program, impowers the users to study 
+! the thermodynamical properties of spins arranged in any complex 
+! lattice network.
+
+! Copyright (C) 2021  Mukesh Kumar Sharma ([email protected])
+
+! 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, see https://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html
+
+	subroutine Monte_Carlo(accept_count)
+
+		use init
+		use omp_lib
+
+		implicit none
+
+		integer :: i, j, k, ith, jth, kth, lth
+		real(8) :: S_vec_previous(5), S_vec_updated(5), S_trial_present(3), &
+			total_eng, eta
+		integer, intent(out) :: accept_count
+
+		accept_count = 0
+		call omp_set_nested(.true.)
+
+		!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i, j, k, ith, jth, kth, lth, &
+		!$OMP& S_vec_previous, S_vec_updated, S_trial_present, total_eng, eta)
+
+		!$OMP DO SCHEDULE(DYNAMIC) COLLAPSE(3)
+		do i = 1, sc(1)
+			do j = 1, sc(2)
+				do k = 1, sc(3)
+
+			ith = bblx(i)
+			jth = bbly(j)   
+			kth = bblz(k)   
+
+			do lth = 1, lattice_per_unit_cell
+
+			! Copy current spin state
+			S_vec_previous(1:5) = ion(1:5, ith, jth, kth, lth)
+
+			! Update spin details
+			call update_spin_details(S_vec_previous, S_vec_updated)
+
+			! Calculate the trial spin
+			S_trial_present(1:3) = S_vec_updated(1:3) - S_vec_previous(1:3)
+
+			! Calculate energy from Hamiltonian
+			call Hamiltonian(ith, jth, kth, lth, S_trial_present, total_eng)
+
+			! Get random number for Metropolis criterion
+			call get_random_num(0d0, 1d0, eta)
+
+			! Metropolis acceptance algorithm
+			if (exp(-beta*total_eng) .gt. eta) then
+			ion(1:5, ith, jth, kth, lth) = S_vec_updated(1:5)
+
+			!$OMP ATOMIC
+			accept_count = accept_count + 1
+			!$OMP END ATOMIC
+
+                        call boundary_condition(ith, jth, kth, lth)
+			end if
+
+			end do
+
+				end do
+			end do
+		end do
+
+		!$OMP END DO
+		!$OMP END PARALLEL
+
+	end subroutine Monte_Carlo
+