diff --git a/src/George_Marsaglia.f90 b/src/George_Marsaglia.f90 new file mode 100644 index 0000000..1797ccf --- /dev/null +++ b/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 (msharma1@ph.iitr.ac.in) + +! 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 diff --git a/src/Hamiltonian.f90 b/src/Hamiltonian.f90 new file mode 100644 index 0000000..4a45962 --- /dev/null +++ b/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 (msharma1@ph.iitr.ac.in) + +! 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 + diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..76ce58b --- /dev/null +++ b/src/Makefile @@ -0,0 +1,3 @@ +ROOTDIR=.. + +include ../Makefile diff --git a/src/Makefile.ether b/src/Makefile.ether new file mode 100644 index 0000000..61a36d1 --- /dev/null +++ b/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 diff --git a/src/Monte_Carlo.f90 b/src/Monte_Carlo.f90 new file mode 100644 index 0000000..9917ccc --- /dev/null +++ b/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 (msharma1@ph.iitr.ac.in) + +! 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 + diff --git a/src/allocate_observables.f90 b/src/allocate_observables.f90 new file mode 100644 index 0000000..d68fd41 --- /dev/null +++ b/src/allocate_observables.f90 @@ -0,0 +1,33 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 allocate_observables + + use init + + implicit none + + allocate(e_mag_avg(sample), e_U_mag(sample), e_chi(sample), & + e_eng_avg(sample), e_U_eng(sample), e_cv(sample), acceptance_ratio(nscan), & + mm_vector(nspecies, 1:3), temp_T(nscan), s_mag_avg_T(nscan), s_chi_T(nscan), & + err_mag_avg_T(nscan), err_chi_T(nscan), s_U_mag_T(nscan), err_U_mag_T(nscan), & + s_eng_avg_T(nscan), s_cv_T(nscan), err_eng_avg_T(nscan), err_cv_T(nscan), & + s_U_eng_T(nscan), err_U_eng_T(nscan), mm_vector_avg_T(nspecies, 1:3, nscan)) + + end subroutine allocate_observables diff --git a/src/boundary_condition.f90 b/src/boundary_condition.f90 new file mode 100644 index 0000000..9ca9912 --- /dev/null +++ b/src/boundary_condition.f90 @@ -0,0 +1,59 @@ +! Code Ether, based on the Monte Carlo technique, can be used to +! study the static and dynamics of spin models applied to any ion geometry. +! Copyright (C) 2021 Mukesh Kumar Sharma (msharma1@ph.iitr.ac.in) + +! 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 boundary_condition(ith, jth, kth, lth) + + use init + + implicit none + + integer, intent(in) :: ith, jth, kth, lth + + ! BOUNDARY X + along_x : if(bc(1).eq.'c')then + if (ith.gt.sc(1)) then + ion(1:5, ith - sc(1), jth, kth, lth) = & + ion(1:5, ith, jth, kth, lth) + elseif(ith.lt.(fromx + nbd_cell_x))then + ion(1:5, ith + sc(1), jth, kth, lth) = & + ion(1:5, ith, jth, kth, lth) + end if + end if along_x + + ! BOUNDARY Y + along_y : if(bc(2).eq.'c')then + if (jth.gt.sc(2)) then + ion(1:5, ith, jth - sc(2), kth, lth) = & + ion(1:5, ith, jth, kth, lth) + elseif(jth.lt.(fromy + nbd_cell_y))then + ion(1:5, ith, jth + sc(2), kth, lth) = & + ion(1:5, ith, jth, kth, lth) + end if + end if along_y + + ! BOUNDARY Z + along_z : if(bc(3).eq.'c')then + if (kth.gt.sc(3)) then + ion(1:5, ith, jth, kth - sc(3), lth) = & + ion(1:5, ith, jth, kth, lth) + elseif(kth.lt.(fromz + nbd_cell_z))then + ion(1:5, ith, jth, kth + sc(3), lth) = & + ion(1:5, ith, jth, kth, lth) + end if + end if along_z + + end subroutine boundary_condition diff --git a/src/ether.f90 b/src/ether.f90 new file mode 100644 index 0000000..ef4befb --- /dev/null +++ b/src/ether.f90 @@ -0,0 +1,307 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 + + program ether + + use init + use omp_lib + use mpi + + implicit none + + ! DATE & TIME + character(8) :: date + character(10) :: time + character(5) :: zone + + integer, dimension(8) :: start, finish + integer :: days, hrs, mins, secs, n_seed, & + stepi, ei, ej, el + integer, dimension(:), allocatable :: seed + + ! Initialize MPI + call MPI_INIT(ierr) + call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr) + call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr) + + call random_seed(size=n_seed) + allocate(seed(n_seed)) + + ! initialize the random number seed = 1992 + seed = 1992 + call random_seed(put=seed) + call date_and_time(VALUES=start) + + if (rank == 0) call system('rm -f *.dat *.xsf') + + call get_tot_species(nspecies, total_ions_per_cell) + + allocate(s(nspecies), sia_vec(1:3, nspecies), & + x(total_ions_per_cell), y(total_ions_per_cell), z(total_ions_per_cell), & + species(0:nspecies), ionn(0:nspecies)) + + call read_input + if (rank == 0) call startup(start) + call read_structure(nspecies, total_ions_per_cell, lp, abc, x, y, z, ionn, species, rank) + call j_values + call get_sia_values + call parameters + call generate_supercell + if (rank == 0) call write_initial_conf + call getting_nbd + if (rank == 0) call write_nbd + call generate_bubble_indices + + ! Total temperature count + nscan = nint((ht - lt) / (tint)) + 1 + + call allocate_observables + + ! For initiating Ground Spin State (GSS) file + if (rank == 0) call write_gss(0) + ! For creating output files + if (rank == 0) call write_output_files(0) + + allocate(addmoreitr(size), num_iterations(size), istart(0:size), iend(0:size), & + tlobs(size), tlspn(size)) + + addmoreitr = 0; istart = 0; iend = 0 + + ! Calculate the reminder (left) values is size divides the nscan, + ! it can be used for assigning this leftovers to other processors. + left = mod(nscan, size) + do el = 1, left + + addmoreitr(el) = 1 ! add more iterations + + end do + + ! Calculate the work distribution: number of iterations per process + interval = nscan / size + do ei = 1, size + + num_iterations(ei) = addmoreitr(ei) + interval + + end do + + ! Calculate the starting and end points of distributed temperature index + do ei = 1, size + + istart(ei) = iend(ei-1) + 1 + iend(ei) = iend(ei-1) + num_iterations(ei) + + end do + + ! Number of observable (like: temp., energy, Cv, Mag., Chi,..., etc.) + total_observables = 14 + + ! Allocate array for each process + ! local observable length + local_olen = total_observables + 3 * nspecies ! total observable + + ! 3 magnetic component * no. of species + local_slen = product(sc) * lattice_per_unit_cell * 4 ! total no. of cell * lattice per unit cel * + ! (3 spin vectors + ID) + + do ei = 1, size + + ! Total length of observables for local array + tlobs(ei) = local_olen*(iend(ei) - istart(ei) + 1) + + ! Total length of spin details for local array + tlspn(ei) = local_slen*(iend(ei) - istart(ei) + 1) + + end do + + allocate(local_obs(tlobs(rank+1)), local_spn(tlspn(rank+1)), si_obs(size), & + si_spn(size)) + + local_obs = 0d0; local_spn = 0d0 + + ! Calculate the starting index for recieving buffer from non-root processors + si_obs(1) = 0 + si_spn(1) = 0 + do ei = 2, size + + ! For observables + si_obs(ei) = sum(tlobs(1:ei-1)) + ! For spin details + si_spn(ei) = sum(tlspn(1:ei-1)) + + end do + + if (rank == 0) then + + !$OMP PARALLEL + num_of_threads = omp_get_num_threads() + !$OMP END PARALLEL + + write(6, *) '' + write(6,'(" Calculations are running on total: ")') + write(6, *) '' + write(6,'(" OpenMP threads =>", i3)') num_of_threads + write(6,'(" MPI processes =>", i3)') size + write(6, *) '' + + write(6, "(' List of temperatures')") + write(6, "(' ~~~~~~~~~~~~~~~~~~~~')") + write(6, *) '' + + do ei = 0, size-1 + + allocate(temp_assigned(num_iterations(ei + 1))) + el = 0 + do ej = istart(ei + 1), iend(ei + 1) + el = el + 1 + temp_assigned(el) = ht - tint * (ej - 1) + end do + write(6, "(' At process ID', i3, ' are:')") ei + write(6, "(' ',5g11.4)") temp_assigned + write(6, *) '' + deallocate(temp_assigned) + + end do + + ! Allocate global array only on the root process (rank 0) + allocate(global_obs(local_olen*nscan), global_spn(local_slen*nscan)) + + else + ! This is necessay because for command line flag -check=all it gives error + ! so we are making these buffer as dummy things for non-root processes + allocate(global_obs(0), global_spn(0)) + end if + + el = 0 + li_obs = 0 ! local index for observables + li_spn = 0 ! local index for spins + temperature: do itemp = istart(rank + 1), iend(rank + 1) + + ! local index for local_* array + el = el + 1 + li_obs = local_olen*(el - 1) + li_spn = local_slen*(el - 1) + + ! Temperature (T) + temp = ht - tint * (itemp - 1) + ! kBT^-1 + beta = 1.0d0 / (kb*temp) + + ! Creating *spK* file + initiate_spin_files = .TRUE. + call write_spins_at_K(itemp) + + ! Setting zero to sample variables + call zeroes('sample') + + acceptance_counting = 0 + sampling: do samplei = 1, sample + + call fresh_spins + call zeroes('eng_mag') + + do stepi = 1, tmcs + + if (ovrr.and..not.Ising.and. & + (mod(real(stepi), real(ovrr_start_interval)) == 0.0)) then + + call overrelaxation + + end if + + call mc(stepi, acceptance_count) + acceptance_counting = acceptance_counting + acceptance_count + + end do + + call evaluate_observables('avg_eng') + call evaluate_observables('avg_mag') + + end do sampling + + call evaluate_observables('observables_per_sample') + + ! At temperature K + + ! Getting magnetic moment vectors + call get_moment_vectors + ! Evaluate all observables + call evaluate_observables('store') + + ! Writing spin states into *spK* file + call write_spins_at_K(itemp) + + ! Writing Job status + write(6, "(' Task for temperature ', g11.4, 'is completed.')") temp + + end do temperature + + ! Gather data from all processes to root + call MPI_GATHERV(local_obs, tlobs(rank+1), MPI_DOUBLE_PRECISION, & + global_obs, tlobs, si_obs, & + MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + + call MPI_GATHERV(local_spn, tlspn(rank+1), MPI_DOUBLE_PRECISION, & + global_spn, tlspn, si_spn, & + MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + + write(6, "(' Process ID', i3, ' is free now, assigned work is completed.')") rank + + ! Root process will finalize observables and write results after all data has been received + if (rank == 0) then + + call evaluate_observables('write') + call graph + + call system('rm -rf data spins') + call system('mkdir spins') + call system('mkdir data') + call system('mv *spK* spins') + call system('mv *conf.xsf *.dat data') + call system('mv graph* data') + + write(6, *) '' + write(6, *) "May the force be with you" + write(6, *) "~ Mukesh Kumar Sharma" + write(6, *) "e-mail@ msharma1@ph.iitr.ac.in" + write(6, *) '' + write(6, *) '' + + close(10004) + close(10005) + + write(6, *) '' + write(6, *) '' + write(6, *) 'SUMMARY :::::::::::::::::::::::::::::::::::::' + write(6, *) '' + write(6, '(" PROGRAM STARTED on date ",i2,"-",i2,"-",i4)') start(3), start(2), start(1) + write(6, "(' at time ',i2,' hrs. ',i2,' min. ',i2,' sec. ')") start(5:7) + call date_and_time(VALUES=finish) + write(6, *) '' + write(6, '(" PROGRAM ENDED on date ",i2,"-",i2,"-",i4)') finish(3), finish(2), finish(1) + write(6, "(' at time ',i2,' hrs. ',i2,' min. ',i2,' sec. ')") finish(5:7) + write(6, *) '' + write(6, *) ':::::::::::::::::::::::::::::::::::::::::::::' + write(6, *) '' + + end if + + ! Finalize MPI + call MPI_FINALIZE(ierr) + + end program ether + diff --git a/src/evaluate_observables.f90 b/src/evaluate_observables.f90 new file mode 100644 index 0000000..16ffdda --- /dev/null +++ b/src/evaluate_observables.f90 @@ -0,0 +1,236 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 evaluate_observables(observable_case) + + use init + + implicit none + + integer :: i, j, k, l, m + real(8) :: mag_value, U_mag, chi, U_eng, cv, sample_1, mag_per_site, & + eng_per_site + + character(len=*), intent(in) :: observable_case + + select case(observable_case) + + case('eng') + + eng_per_site = eng/(2*total_lattice_sites) + eng_avg = eng_avg + eng_per_site + eng2_avg = eng2_avg + eng_per_site**2.0d0 + eng4_avg = eng4_avg + eng_per_site**4.0d0 + + case('mag') + + mag_value = sqrt(dot_product(net_mag, net_mag)) + mag_per_site = mag_value/total_lattice_sites + mag_avg = mag_avg + mag_per_site + mag2_avg = mag2_avg + mag_per_site**2.0d0 + mag4_avg = mag4_avg + mag_per_site**4.0d0 + + case('avg_mag') + + mag_avg = mag_avg/total_calculations + mag2_avg = mag2_avg/total_calculations + mag4_avg = mag4_avg/total_calculations + s_mag_avg = mag_avg + s_mag_avg + e_mag_avg(samplei) = mag_avg + + U_mag = 1.0d0 - ((1.0d0/3.0d0)*(mag4_avg/(mag2_avg**2.0d0))) + s_U_mag = U_mag + s_U_mag + e_U_mag(samplei) = U_mag + chi = beta*(mag2_avg - mag_avg**2.0d0)*total_lattice_sites + s_chi = chi + s_chi + e_chi(samplei) = chi + + case('avg_eng') + + eng_avg = eng_avg/total_calculations + eng2_avg = eng2_avg/total_calculations + eng4_avg = eng4_avg/total_calculations + + s_eng_avg = eng_avg + s_eng_avg + e_eng_avg(samplei) = eng_avg + + U_eng = 1.0d0 - ((1.0d0/3.0d0)*(eng4_avg/(eng2_avg**2))) + s_U_eng = U_eng + s_U_eng + e_U_eng(samplei) = U_eng + + cv = (beta**2.0d0)*(eng2_avg - eng_avg**2.0d0)*total_lattice_sites + s_cv = cv + s_cv + e_cv(samplei) = cv + + case('observables_per_sample') + + sample_1 = 1.0d0/(sample*(sample - 1.0d0)) + + !MAGNETIZATION PER SAMPLE + s_mag_avg = s_mag_avg/sample + s_U_mag = s_U_mag/sample + s_chi = s_chi/sample + + ! ERRORS + err_U_mag = 0.0; err_chi = 0.0; err_mag_avg = 0.0 + + do i = 1, sample + + err_U_mag = err_U_mag + (e_U_mag(i) - s_U_mag)**2 + err_chi = err_chi + (e_chi(i) - s_chi)**2 + err_mag_avg = err_mag_avg + (e_mag_avg(i) - s_mag_avg)**2 + + end do + + err_U_mag = sqrt(sample_1*err_U_mag) + err_chi = sqrt(sample_1*err_chi) + err_mag_avg = sqrt(sample_1*err_mag_avg) + + !ENERGY PER SAMPLE + s_eng_avg = s_eng_avg/sample + s_U_eng = s_U_eng/sample + s_cv = s_cv/sample + + ! ERRORS + err_U_eng = 0.0; err_cv = 0.0; err_eng_avg = 0.0 + + do i = 1, sample + + err_U_eng = err_U_eng + (e_U_eng(i) - s_U_eng)**2 + err_cv = err_cv + (e_cv(i) - s_cv)**2 + err_eng_avg = err_eng_avg + (e_eng_avg(i) - s_eng_avg)**2 + + end do + + err_U_eng = sqrt(sample_1*err_U_eng) + err_cv = sqrt(sample_1*err_cv) + err_eng_avg = sqrt(sample_1*err_eng_avg) + + case('store') + + ! ENERGY + !temp_T(itemp) = temp; s_mag_avg_T(itemp) = s_mag_avg; s_chi_T(itemp) = s_chi + !err_mag_avg_T(itemp) = err_mag_avg; err_chi_T(itemp) = err_chi; s_U_mag_T(itemp) = s_U_mag + !err_U_mag_T(itemp) = err_U_mag + + local_obs(1 + li_obs) = temp + local_obs(2 + li_obs) = s_mag_avg + local_obs(3 + li_obs) = s_chi + local_obs(4 + li_obs) = err_mag_avg + local_obs(5 + li_obs) = err_chi + local_obs(6 + li_obs) = s_U_mag + local_obs(7 + li_obs) = err_U_mag + + ! MAGNETIC +! s_eng_avg_T(itemp) = s_eng_avg; s_cv_T(itemp) = s_cv; err_eng_avg_T(itemp) = err_eng_avg +! err_cv_T(itemp) = err_cv; s_U_eng_T(itemp) = s_U_eng; err_U_eng_T(itemp) = err_U_eng + + local_obs(8 + li_obs) = s_eng_avg + local_obs(9 + li_obs) = s_cv + local_obs(10 + li_obs) = err_eng_avg + local_obs(11 + li_obs) = err_cv + local_obs(12 + li_obs) = s_U_eng + local_obs(13 + li_obs) = err_U_eng + + ! ACCEPTANCE RATIOS + local_obs(14 + li_obs) = acceptance_counting/(tmcs*total_lattice_sites*sample*1.0)*100 + + j = 0 + do i = 1, nspecies + + j = j + 1 + local_obs(total_observables + (3*(i-1) + 1) + li_obs) = mm_vector(j, 1) + local_obs(total_observables + (3*(i-1) + 1) + 1 + li_obs) = mm_vector(j, 2) + local_obs(total_observables + (3*(i-1) + 1) + 2 + li_obs) = mm_vector(j, 3) + + end do + + m = 0 + do i = fromx, tox + do j = fromy, toy + do k = fromz, toz + do l = 1, lattice_per_unit_cell + + local_spn(m + 1 + li_spn) = ion(1, i, j, k, l)*s(int(ion(4, i, j, k, l))) + local_spn(m + 2 + li_spn) = ion(2, i, j, k, l)*s(int(ion(4, i, j, k, l))) + local_spn(m + 3 + li_spn) = ion(3, i, j, k, l)*s(int(ion(4, i, j, k, l))) + local_spn(m + 4 + li_spn) = ion(4, i, j, k, l) + m = m + 4 + + end do + end do + end do + end do + + case('write') + + lobs = 0 ! local observable lenght + lspn = 0 ! local spin length + do itemp = 1, nscan + + lobs = (itemp -1)*local_olen + lspn = (itemp -1)*local_slen + + temp_T(itemp) = global_obs(1 + lobs) + s_mag_avg_T(itemp) = global_obs(2 + lobs) + s_chi_T(itemp) = global_obs(3 + lobs) + err_mag_avg_T(itemp) = global_obs(4 + lobs) + err_chi_T(itemp) = global_obs(5 + lobs) + s_U_mag_T(itemp) = global_obs(6 + lobs) + err_U_mag_T(itemp) = global_obs(7 + lobs) + + s_eng_avg_T(itemp) = global_obs(8 + lobs) + s_cv_T(itemp) = global_obs(9 + lobs) + err_eng_avg_T(itemp) = global_obs(10 + lobs) + err_cv_T(itemp) = global_obs(11 + lobs) + s_U_eng_T(itemp) = global_obs(12 + lobs) + err_U_eng_T(itemp) = global_obs(13 + lobs) + + acceptance_ratio(itemp) = global_obs(14 + lobs) + + j = 0 + + do i = 1, nspecies + + j = j + 1 + mm_vector_avg_T(j, 1, itemp) = global_obs(total_observables + (3*(i-1) + 1) + lobs) + mm_vector_avg_T(j, 2, itemp) = global_obs(total_observables + (3*(i-1) + 1) + 1 + lobs) + mm_vector_avg_T(j, 3, itemp) = global_obs(total_observables + (3*(i-1) + 1) + 2 + lobs) + + + end do + + ! Writing output files + call write_output_files(itemp) + ! Writing Ground Spin States (GSS) file + call write_gss(itemp) + + end do + + case default + + write(6, *) '' + write(6, "(' ==> Found unknown case tag :',A8 )") observable_case + write(6, *) " in 'evaluate_observables' subroutine" + write(6, *) ' STOPPING now' + stop + + end select + + end subroutine evaluate_observables diff --git a/src/fresh_spins.f90 b/src/fresh_spins.f90 new file mode 100644 index 0000000..e49979c --- /dev/null +++ b/src/fresh_spins.f90 @@ -0,0 +1,46 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 fresh_spins + + use init + + implicit none + + integer :: i, j, k, l + real(8) :: S_vec_previous(5), S_vec_updated(5) + + S_vec_previous = 0; S_vec_updated = 0 + + do i = fromx, tox + do j = fromy, toy + do k = fromz, toz + do l = 1, lattice_per_unit_cell + + S_vec_previous = ion(1:5, i, j, k, l) + + call update_spin_details(S_vec_previous, S_vec_updated) + ion(1:5, i, j, k, l) = S_vec_updated + + end do + end do + end do + end do + + end subroutine fresh_spins diff --git a/src/generate_bubble_indices.f90 b/src/generate_bubble_indices.f90 new file mode 100644 index 0000000..dc69e99 --- /dev/null +++ b/src/generate_bubble_indices.f90 @@ -0,0 +1,61 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 generate_bubble_indices + + use init + + implicit none + + allocate(bblx(sc(1)), bbly(sc(2)), bblz(sc(3))) + + bblx = 0; bbly = 0; bblz = 0 + + call get_index(bblx, nbd_cell_x + 1, fromx, tox) + call get_index(bbly, nbd_cell_y + 1, fromy, toy) + call get_index(bblz, nbd_cell_z + 1, fromz, toz) + + contains + + subroutine get_index(indices, stepi, from, to) + + implicit none + + integer :: i, j, k, extension + integer, intent(in) :: stepi, from, to + integer, allocatable, intent(out) :: indices(:) + + extension = to - from + 1 + allocate(indices(extension)) + j = from; k = 0 + do i = 1, extension + + indices(i) = j + stepi*k + k = k + 1 + if(indices(i).gt.to) then + j = j + 1; k = 0 + indices(i) = j + k = k + 1 + end if + + end do + + end subroutine get_index + + end subroutine generate_bubble_indices diff --git a/src/generate_supercell.f90 b/src/generate_supercell.f90 new file mode 100644 index 0000000..0d74b76 --- /dev/null +++ b/src/generate_supercell.f90 @@ -0,0 +1,168 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 generate_supercell + + use init + + implicit none + + integer :: i, j, k, l, s_count + + real(8), allocatable :: atom(:, :), ar(:, :) + real(8) :: sx, sy, sz, rnum + logical :: file_found + + ! SAVING DETAILS OF ATOMS PRESENT IN STRUCTURE FILE + allocate(atom(total_ions_per_cell, 0:8)) + atom = 0 + atomic_details: do k = 1, nspecies + do i = int(sum(ionn(0:k-1)))+1, int(sum(ionn(0:k))) + + atom(i,0) = i ! atom no. + if(Ising) then + atom(i, 3) = (-1)**(i) + else + if (angle.eqv..FALSE.) then + call George_Marsaglia(sx, sy, sz) + atom(i, 3) = sz ! S(z) + atom(i, 1) = sx ! S(x) + atom(i, 2) = sy ! S(y) + atom(i, 5) = 0 ! phi value + else + + call get_random_num(-1d0, +1d0, rnum) + atom(i, 3) = rnum ! S(z) + call get_random_num(0.0d0, 2d0*pi, rnum) + atom(i, 1) = sqrt(1d0-(atom(i, 3)**2))*cos(rnum) ! S(x) + atom(i, 2) = sqrt(1d0-(atom(i, 3)**2))*sin(rnum) ! S(y) + atom(i, 5) = rnum ! phi value + end if + endif + atom(i, 4) = k ! tag for element + atom(i, 6) = x(i) ! x + atom(i, 7) = y(i) ! y + atom(i, 8) = z(i) ! z + + enddo + enddo atomic_details + + ! FILTRATION OF SPECIES FOR CALC. + s_count = 0 + do i = 1, n_speci_incl + do j =1 , nspecies + if (specicies_to_include(i).eq.species(j)) then + s_count = s_count + ionn(j) + end if + end do + end do + + allocate(ar(s_count, 0:8), stgg_ion(s_count)) + ar = 0 + + ! STAGGERED INFO + + if (staggered) then !getting staggered values of each ion + inquire(file='staggered', exist=file_found) + + if(file_found) then + if (rank == 0) write(6, *) '' + if (rank == 0) write(6, *) ' :::::::::::::::::::::::::: & + STG list ::::::::::::::::::::::::' + open(unit=2, file='staggered', status='old', action='read') + do i = 1, s_count + read(2, *) stgg_ion(i) + end do + close(2) + if (rank == 0) write(6, "(/,' > ',8f7.1)") stgg_ion + if (rank == 0) write(6, *) '' + else + if (rank == 0) then + write(6, *) '' + write(6, *) ' Logic for staggered magnetiaztion is .TRUE.' + write(6, *) " but 'staggered' file is not provided" + write(6, *) ' ~~~~~~~~~' + write(6, *) '' + write(6, *) " file formate for 'staggered' is as follows" + write(6, *) '' + write(6, *) ' c1 for atom1' + write(6, *) ' c2 for atom1' + write(6, *) ' c3 for atom3' + write(6, *) ' . .' + write(6, *) ' . .' + write(6, *) ' . .' + write(6, *) ' Note: Sequence of c[i] for ith atom should' + write(6, *) ' follows the sequnce as according to' + write(6, *) ' provided structure file' + write(6, *) '' + write(6, *) ' STOPPING now' + write(6, *) '' + stop + end if + end if + end if + + ! Getting total lattice per unit cell, which will be used for MC calculations + lattice_per_unit_cell = 0 + inclusion : do i = 1, n_speci_incl + do j =1 , nspecies + if (specicies_to_include(i).eq.species(j)) then + + do k = sum(ionn(0:j-1))+1, sum(ionn(0:j)) + lattice_per_unit_cell = lattice_per_unit_cell +1 ! Counting included ions + do l = 0, 8 + ar(lattice_per_unit_cell, l) = atom(k, l) ! 'k' will be used for identifying + end do ! these element in further last moment + enddo + endif + enddo + enddo inclusion + deallocate(atom) + + ! EXPANDING INTO SUPERCELL + allocate(ion(0:8, sc(1) + 2*nbd_cell_x, sc(2) + 2*nbd_cell_y, & + sc(3) + 2*nbd_cell_z, lattice_per_unit_cell)) + s_count = 0; ion = 0 + do i = 1, sc(1) + 2*nbd_cell_x + do j = 1, sc(2) + 2*nbd_cell_y + do k = 1, sc(3) + 2*nbd_cell_z + do l = 1, lattice_per_unit_cell + + s_count = s_count +1 ! Counting total no. of ions in supercell + ion(0, i, j, k, l) = s_count ! ID + ion(1:5, i, j, k, l) = ar(l, 1:5) ! Sx, Sy, Sz, species no. & phi + ion(6:8, i, j, k, l) = abc(1, 1:3)*(i - 1) & ! co-ordinates + + abc(2, 1:3)*(j - 1) + abc(3, 1:3)*(k - 1) & + + ar(l, 6:8) + + end do + end do + end do + end do + + deallocate(ar) + + ! total no. of lattice points + total_ions = s_count ! in super cell + total_lattice_sites = product(sc)*lattice_per_unit_cell ! in super cell simulation box + + if (rank == 0) write(6, *) '==> Supercells are generated' + + end subroutine generate_supercell + diff --git a/src/get_moment_vectors.f90 b/src/get_moment_vectors.f90 new file mode 100644 index 0000000..33d45c5 --- /dev/null +++ b/src/get_moment_vectors.f90 @@ -0,0 +1,45 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 get_moment_vectors + + use init + + implicit none + + integer :: i, j, k, l + + mm_vector = 0 + + do k = fromz, toz + do j = fromy, toy + do i = fromx, tox + do l = 1, lattice_per_unit_cell + + ! moment vectors for each magentic ion + mm_vector(int(ion(4, i, j, k, l)), 1:3) = & + mm_vector(int(ion(4, i, j, k, l)), 1:3) & + + ion(1:3, i, j, k, l)*s(int(ion(4, i, j, k, l))) + + end do + end do + end do + end do + + end subroutine get_moment_vectors diff --git a/src/get_ovrr_vec.f90 b/src/get_ovrr_vec.f90 new file mode 100644 index 0000000..b4fbf72 --- /dev/null +++ b/src/get_ovrr_vec.f90 @@ -0,0 +1,133 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 get_ovrr_vec(io, jo, ko, lo, ovrr_vec) + + use init + + implicit none + + integer, intent(in) :: io, jo, ko, lo + integer :: ionID, ith_bond, ith_nbr, cell + + real(8) :: Si_ref(3), sia_value(3) + real(8), intent(inout) :: ovrr_vec(3) + + ovrr_vec = 0 + + ! Central ion's species ID + ionID = int(ion(4, io, jo, ko, lo)) + + ! Central ion's spin vectors + Si_ref = ion(1:3, io, jo, ko, lo) + + call JSiSj_ovrr(io, jo, ko, lo, ionID, ovrr_vec) + + if(Zeeman) call gmbSH_ovrr(g_factor, mb, H, ovrr_vec) + + if(single_ion_anisotropy) then + sia_value(1:3) = sia_vec(1:3, ionID) + call siaS2_ovrr(Si_ref, sia_value, ovrr_vec) + end if + + contains + + ! JSiSj + subroutine JSiSj_ovrr(i, j, k, l, central_ion_ID, ovrr_vec) + + implicit none + + integer, intent(in) :: i, j, k, l, central_ion_ID + integer :: total_nbr, posx, posy, posz, ion_ID + + real(8) :: Sj(3), Jij(3), JSj(3) + real(8), intent(inout) :: ovrr_vec(3) + + ion_ID = int(ion(0, i, j, k, l)) + + ! Over distinct bonds + do ith_bond = 1, no_of_nbd + + ! For ith_bond bond total connecting neighbours to the central ion [ID = ion(0, i, j, k, l)] is + total_nbr = nn(ith_bond, ion_ID, 0, 0) + + ! no. of similar nbd for ith distinct bond (ith_bond) + do ith_nbr = 1, total_nbr + + ! nbr's cell position + posx = nn(ith_bond, ion_ID, ith_nbr, 1) + posy = nn(ith_bond, ion_ID, ith_nbr, 2) + posz = nn(ith_bond, ion_ID, ith_nbr, 3) + + ! nbr's ID in the cell + cell = nn(ith_bond, ion_ID, ith_nbr, 4) + + if(ion(4, posx, posy, posz, cell).eq.0) then + go to 4 + end if + + Sj(1:3) = ion(1:3, posx, posy, posz, cell) + + !Jij term + Jij = j_exc(ith_bond, central_ion_ID, int(ion(4, posx, posy, posz, cell)), 1:3) + + !Jij.Sj term + JSj = Jij*Sj + + ovrr_vec = ovrr_vec + JSj + +4 continue + end do + end do + + end subroutine JSiSj_ovrr + + !gmbSH + subroutine gmbSH_ovrr(g, mu, H, ovrr_vec) + + implicit none + + real(8) :: ref_vec(3) + real(8), intent(in) :: g, mu, H(3) + real(8), intent(inout) :: ovrr_vec(3) + + ! OVRR vec. due to magnetic field + ref_vec = -(g*mu*H) + + ovrr_vec = ovrr_vec + ref_vec + + end subroutine gmbSH_ovrr + + !SINGLE ION ANISOTRPY (SIA) + subroutine siaS2_ovrr(Si, sia_val, ovrr_vec) + + implicit none + + real(8) :: ref_vec(3) + real(8), intent(in) :: Si(3), sia_val(3) + real(8), intent(inout) :: ovrr_vec(3) + + ! OVVR vec. due to single ion anisiatropy + ref_vec = Si*sia_val + + ovrr_vec = ovrr_vec + ref_vec + + end subroutine siaS2_ovrr + + end subroutine get_ovrr_vec diff --git a/src/get_random_num.f90 b/src/get_random_num.f90 new file mode 100644 index 0000000..6fe9e07 --- /dev/null +++ b/src/get_random_num.f90 @@ -0,0 +1,32 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 get_random_num(from, to, rn) + + use init ! for using seed variable + + implicit none + + real(8), intent(in) :: from, to + real(8), intent(out) :: rn + + call random_number(rn) + rn = (to-from)*rn + from + + end subroutine diff --git a/src/get_sia_values.f90 b/src/get_sia_values.f90 new file mode 100644 index 0000000..92aa23d --- /dev/null +++ b/src/get_sia_values.f90 @@ -0,0 +1,77 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 get_sia_values + + use init + + implicit none + + integer :: i, j, SIA_ID + + character(len=3) :: atom1, ab, out + + logical :: file_found + + if(single_ion_anisotropy.and..not.Ising) then + inquire(file='single_ion_anisotropy', exist=file_found) + if(file_found) then + open(10011, file='single_ion_anisotropy', status='old', action='read') + if (rank == 0) write(6, '(3X,":::::::::: SIA list (meV) ::::::::::")') + if (rank == 0) write(6, *) '' + do i = 1, nspecies + read(10011, *) atom1, sia(1:3) ! Single Ion Anisotropy (SIA) vector in unit of meV + sia_vec(1:3, i) = sia(1:3) + ab = atom1 + call lu(ab(1:1), out, 'U' ) + atom1(1:1) = out + call lu(ab(2:2), out, 'L' ) + atom1(2:2) = out + + if (rank == 0) write(6, "(4X,A4,3f10.3)") atom1, sia_vec(1:3, i) + + do j = 1, nspecies + if (atom1.eq.species(j)) then + SIA_ID = j + end if + end do + + sia_vec(1:3, SIA_ID) = sia_vec(1:3, SIA_ID) * s(SIA_ID)**2 + end do + write(6, *) '' + close(10011) + else + if (rank == 0) then + write(6, *) " SIA is .TRUE. but file 'single_ion_anisotropy'" + write(6, *) " is not found so, stopping now" + write(6, *) '' + write(6, *) " 'single_ion_anisotropy' file formate are as follows" + write(6, *) " ~~~~~~~~~~~~~~~~~~~~~" + write(6, *) " species1_name 1_SIAx 1_SIAy 1_SIAz" + write(6, *) " species2_name 2_SIAx 2_SIAy 2_SIAz" + write(6, *) " ." + write(6, *) " ." + write(6, *) " ." + write(6, *) '' + stop + end if + end if + end if + + end subroutine get_sia_values diff --git a/src/get_tot_energy.f90 b/src/get_tot_energy.f90 new file mode 100644 index 0000000..a75f115 --- /dev/null +++ b/src/get_tot_energy.f90 @@ -0,0 +1,58 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 get_tot_energy(toten) + + use init + use omp_lib + + implicit none + + integer :: i, j, k, l + real(8), intent(out) :: toten + real(8) :: partial_toten, Spin_vec(3) + + toten = 0.0 + call omp_set_nested(.true.) + + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i, j, k, l, Spin_vec, partial_toten) REDUCTION(+:toten) + !$OMP DO SCHEDULE(DYNAMIC) COLLAPSE(3) + do k = fromz, toz + do j = fromy, toy + do i = fromx, tox + + do l = 1, lattice_per_unit_cell + ! Get spin vector from the ion array + Spin_vec(1:3) = ion(1:3, i, j, k, l) + + ! Calculate partial energy for this spin + call Hamiltonian(i, j, k, l, Spin_vec, partial_toten) + + ! Accumulate the total energy + toten = toten + partial_toten + end do + + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + + end subroutine get_tot_energy + diff --git a/src/get_tot_magnetisation.f90 b/src/get_tot_magnetisation.f90 new file mode 100644 index 0000000..7d9b7de --- /dev/null +++ b/src/get_tot_magnetisation.f90 @@ -0,0 +1,69 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 get_tot_magnetisation(magnetisation) + + use init + use omp_lib + + implicit none + + integer :: i, j, k, l, stg + real(8), intent(out) :: magnetisation(3) + real(8) :: temp_magnetisation(3) + + magnetisation = 0.0 + call omp_set_nested(.true.) + + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i, j, k, l, stg, temp_magnetisation) & + !$OMP& REDUCTION(+:magnetisation) ! Reduction on all elements of the magnetisation array + + temp_magnetisation = 0.0 + + !$OMP DO SCHEDULE(DYNAMIC) COLLAPSE(3) + do k = fromz, toz + do j = fromy, toy + do i = fromx, tox + + do l = 1, lattice_per_unit_cell + + ! Get the stg value depending on the 'staggered' flag + if (staggered .eqv. .TRUE.) then + stg = stgg_ion(l) + else + stg = 1 + end if + + temp_magnetisation(1:3) = temp_magnetisation(1:3) + & + stg*ion(1:3, i, j, k, l)*s(int(ion(4, i, j, k, l))) + + end do + + end do + end do + end do + !$OMP END DO + + ! Accumulating net magnetisation + magnetisation(1:3) = magnetisation(1:3) + temp_magnetisation(1:3) + + !$OMP END PARALLEL ! End of the parallel region + +end subroutine get_tot_magnetisation + diff --git a/src/get_tot_species.f90 b/src/get_tot_species.f90 new file mode 100644 index 0000000..75283d7 --- /dev/null +++ b/src/get_tot_species.f90 @@ -0,0 +1,76 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 get_tot_species(n, total_ions_in_cell) + + implicit none + + integer :: i, j, k, jj + integer, intent(out) :: n, total_ions_in_cell + integer, allocatable :: totion(:) + + character(len=200) :: titel + + logical :: file_found + + inquire(file='structure.vasp', exist=file_found) + if(.not.file_found) then + write(6, *) "==> 'structure.vasp' is not present" + write(6, *) " STOPPING now" + write(6, *) "" + stop + end if + + open(unit=0, file='structure.vasp', status='old', action='read') + + n = 0 + + j = 0; k = 0 + do jj = 1, 6 + read(0, *) ! skipping 6 lines + end do + + read(0, '(a)') titel + + i = 1 +100 continue + check: if (titel(i:i).ne.' ') then + k = i+1 +102 continue + if(titel(k:k) .ne.' ')then + k = k+1 + goto 102 + else + j = j + 1 + end if + i = k + 1 + else + i = i + 1 + end if check + + if (i.lt.(len_trim(titel)+1)) goto 100 + n = j !no. of species + + allocate(totion(n)) + read(titel, *) totion ! total ions + total_ions_in_cell = sum(totion) + + close(0) + + end subroutine get_tot_species diff --git a/src/getting_nbd.f90 b/src/getting_nbd.f90 new file mode 100644 index 0000000..a07df6f --- /dev/null +++ b/src/getting_nbd.f90 @@ -0,0 +1,122 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 getting_nbd + + use init + use omp_lib + + implicit none + + integer :: search_nbd_x(2), search_nbd_y(2), search_nbd_z(2), i, j, k, l, & + ii, jj, kk, ll, nbd_count, nbd_dis_max, central_ion_ID, nbd_ID, ith_bond + + real(8) :: distance, origin(3), nbd_pos(3) + + allocate(nn(no_of_nbd, total_ions, 0:20, 0:4)) + + if (rank == 0) write(6, *) '==> Sensing neighbourhoods' + nn = 0 + +!$OMP PARALLEL DEFAULT(PRIVATE) SHARED(sc, lattice_per_unit_cell, no_of_nbd, ion, nn, nbd_dis, & +!$OMP& nbd_cell_x, nbd_cell_y, nbd_cell_z, nbd_finding_criteria) + !$OMP DO SCHEDULE(dynamic) + central_atom : do k = 1, sc(3) + 2*nbd_cell_z + call nbd_finding_limit(sc(3), nbd_cell_z, k, search_nbd_z) + + do j = 1, sc(2) + 2*nbd_cell_y + call nbd_finding_limit(sc(2), nbd_cell_y, j, search_nbd_y) + + do i = 1, sc(1) + 2*nbd_cell_x + call nbd_finding_limit(sc(1), nbd_cell_x, i, search_nbd_x) + + do l = 1, lattice_per_unit_cell + + central_ion_ID = int(ion(0, i, j, k, l)) + origin = ion(6:8, i, j, k, l) + + ! no. of distinct bond + nbd : do ith_bond = 1, no_of_nbd + nbd_count = 0.0 + + do kk = search_nbd_z(1), search_nbd_z(2) + do jj = search_nbd_y(1), search_nbd_y(2) + do ii = search_nbd_x(1), search_nbd_x(2) + do ll = 1, lattice_per_unit_cell + + nbd_pos = ion(6:8, ii, jj, kk, ll) + nbd_ID = int(ion(0, ii, jj, kk, ll)) + + if(central_ion_ID.ne.nbd_ID)then + + distance = sqrt(dot_product(origin - nbd_pos, origin - nbd_pos )) + + if(abs(nbd_dis(ith_bond) - (distance)).le.nbd_finding_criteria)then + + nbd_count = nbd_count +1 + ! Storing nbd's ID + nn(ith_bond, central_ion_ID, nbd_count, 0) = nbd_ID + + ! Storing nbd's indices + nn(ith_bond, central_ion_ID, nbd_count, 1) = ii + nn(ith_bond, central_ion_ID, nbd_count, 2) = jj + nn(ith_bond, central_ion_ID, nbd_count, 3) = kk + + ! Storing nbd's atom count + nn(ith_bond, central_ion_ID, nbd_count, 4) = ll + + end if + end if + + end do + end do + end do + end do + + ! Storing total no. of similar + ! nbds on particular ion for mth distinct bond + nn(ith_bond, central_ion_ID, 0, 0) = nbd_count + + end do nbd + end do + end do + end do + end do central_atom + !$OMP END DO +!$OMP END PARALLEL + if (rank == 0) write(6, *) '==> Neighbourhood informations are created now' + + contains + + subroutine nbd_finding_limit(total_extension, nbd_upto, cell_pos, to_nbd_cell) + + implicit none + + integer, intent (in) :: total_extension, nbd_upto, cell_pos + integer, intent (out) :: to_nbd_cell(2) + + to_nbd_cell(1) = cell_pos - nbd_upto + to_nbd_cell(2) = cell_pos + nbd_upto + if(cell_pos .le. nbd_upto) to_nbd_cell(1) = 1 + if(cell_pos .gt. (total_extension + 2*nbd_upto - nbd_upto)) & + to_nbd_cell(2) = total_extension + 2*nbd_upto + + end subroutine nbd_finding_limit + + end subroutine getting_nbd diff --git a/src/graph.f90 b/src/graph.f90 new file mode 100644 index 0000000..e65c285 --- /dev/null +++ b/src/graph.f90 @@ -0,0 +1,103 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 graph + + use init + + implicit none + + if(para)then + lbl = 'k_{B}T/J' + else + lbl = "K" + end if + + open(10008, file='graph.sh', status='unknown') + write(10008,*) '# set labels, lables, and xtics/ytics accordingly' + write(10008,*) 'set nokey' + write(10008,*) 'set terminal png size 1100, 900 font "Times-New-Roman,18"' + write(10008,*) "set output 'results_Ether.png'" + write(10008,*) 'set multiplot' + write(10008,*) "set size squar 0.5,0.5" + write(10008,*) 'set nokey' + write(10008,*) 'set format y "%g"' + write(10008, 303) lt - tint, abs(ht - lt + tint)/5., ht + write(10008,*) 'set mxtics 10' + write(10008,*) 'set origin 0.0,0.0' + write(10008,*) 'set grid' + if(para)then + write(10008, 301) trim(lbl) + else + write(10008, 302) trim(lbl) + end if + if(staggered) then + write(10008,*) 'set ylabel "Staggered Magnetization (M)"' + else + write(10008,*) 'set ylabel "Magnetization (M)"' + end if + + write(10008,*) "plot 'magnetization.dat' u 1:2:4 with yerrorbars & + lt 6 notitle" + write(10008,*) 'set origin 0.5, 0.0' + if(para)then + write(10008, 301) lbl + else + write(10008, 302) lbl + end if + + if(staggered) then + write(10008,*) 'set ylabel "Staggered Susceptibility (χ)"' + else + write(10008,*) 'set ylabel "Susceptibility (χ)"' + end if + + write(10008,*) "plot 'magnetization.dat' u 1:3:5 with yerrorbars & + lt 6 notitle" + write(10008,*) 'set origin 0.0,0.5' + if(para)then + write(10008, 301) lbl + else + write(10008, 302) lbl + end if + write(10008,*) 'set ylabel "Energy (meV)"' + write(10008,*) "plot 'energy.dat' u 1:($2*1000):4 with yerrorbars & + lt 6 notitle" + write(10008,*) 'set origin 0.5,0.5' + if(para)then + write(10008, 301) lbl + else + write(10008, 302) lbl + end if + if(para)then + write(10008, *) 'set ylabel "Specific Heat (C_{V})"' + else + write(10008, *) 'set ylabel "Specific Heat & + (C_{V}k_{B}^{-1})"' + end if + write(10008,*) "plot 'energy.dat' u 1:3:5 with yerrorbars & + lt 6 notitle" + write(10008,*) 'unset multiplot' + close (10008) + +301 format(' set xlabel "Temperature (',A8,')"') +302 format(' set xlabel "Temperature (',A1,')"') +303 format (" set xtics ", f11.5,",", f11.5,",", f11.5) + + end subroutine graph diff --git a/src/init.f90 b/src/init.f90 new file mode 100644 index 0000000..e005656 --- /dev/null +++ b/src/init.f90 @@ -0,0 +1,78 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 + + module init + + use mpi + + implicit none + + integer, public, save :: tmcs, tmcs_eq, n_speci_incl, sc(3), sample, samplei, & + overrelaxed, ovrr_start_interval, optbeta, nspecies, total_ions_per_cell, & + no_of_nbd, j_ID(2), similar_bonds, num_of_threads, total_lattice_sites, & + lattice_per_unit_cell, nscan, nbd_cell_x, nbd_cell_y, nbd_cell_z, & + total_ions, fromx, fromy, fromz, tox, toy, toz, total_calculations, itemp, & + spin_file_ID, gss_ID, acceptance_counting, acceptance_count + + real(8), public, parameter :: pi = 4*atan(1d0) + + real(8), public, save :: temp, ht, lt, tint, dphi, to_cal, g_factor, para_value, & + exchange_interval, lp(3), abc(3, 3), sia(1:3), kb = 8.617333262d-5, & + h(3), anisotropy(3), j_value, mb, beta_critria, & + nbd_finding_criteria = 0.0001d0, convert_to_rad = pi/180.0d0, beta, & + eng, eng_avg, eng2_avg, eng4_avg, mag_avg, mag2_avg, mag4_avg, & + s_eng_avg, s_eng2_avg, e_eng2_avg, s_U_eng, s_cv, & + s_mag_avg, s_mag2_avg, e_mag2_avg, s_U_mag, s_chi, & + net_mag(3), err_U_mag, err_chi, err_mag_avg, err_U_eng, err_cv, & + err_eng_avg + + integer, allocatable, public, save :: ionn(:), nn(:,:,:,:), bblx(:), bbly(:), bblz(:) + + real(8), allocatable, public, save :: x(:), y(:), z(:), s(:), nbd_dis(:), & + j_exc(:,:,:,:), sia_vec(:, :), stgg(:), stgg_ion(:), ion(:, :, :, :, :), & + e_mag_avg(:), e_U_mag(:), e_chi(:), e_eng_avg(:), e_U_eng(:), e_cv(:), & + mm_vector(:, :), & + temp_T(:), s_mag_avg_T(:), s_chi_T(:), err_mag_avg_T(:), err_chi_T(:), & + s_U_mag_T(:), err_U_mag_T(:), s_eng_avg_T(:), s_cv_T(:), err_eng_avg_T(:), & + err_cv_T(:), s_U_eng_T(:), err_U_eng_T(:), mm_vector_avg_T(:, :, :), & + acceptance_ratio(:) + + character, public, save :: model, bc(3) + character*30, public, save :: title, coordinate, filename, lbl + character(len=2), allocatable, public, save :: specicies_to_include(:), species(:) + character(len=20), dimension(50), public, save :: m_head, e_head + + logical, public, save :: staggered, angle, Zeeman, single_ion_anisotropy, para, ovrr, & + EXalgo, temp_ex, beta_file, initiate_spin_files + +! For MPI's + integer :: rank, size, ierr, interval, left, tag, local_olen, local_slen, & + total_observables, status(MPI_STATUS_SIZE), request_obs, request_spn, & + status_obs, status_spn, li_obs, li_spn, lobs, lspn + + integer, allocatable :: addmoreitr(:), num_iterations(:), & + istart(:), iend(:), tlobs(:), tlspn(:), si_obs(:), & + si_spn(:) + + real(8), allocatable, public, save :: local_obs(:), local_spn(:), & + global_obs(:), global_spn(:), temp_assigned(:) + + logical, public, save :: completed, Ising, Heisenberg + + end module init diff --git a/src/j_values.f90 b/src/j_values.f90 new file mode 100644 index 0000000..365928f --- /dev/null +++ b/src/j_values.f90 @@ -0,0 +1,101 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 j_values + + use init + + implicit none + + integer :: i, j, k, ii, nbd_dis_max + + character(len=3) :: atom(2), atom1, atom2, ab, out + + open(10002, file='j_exchange', status='old', action='read') + read(10002, *) no_of_nbd ! no. of diff. bond length w.r.t. any central ion + allocate(nbd_dis(no_of_nbd)) ! nbd distances + allocate(j_exc(no_of_nbd, nspecies, nspecies, 1:3)) ! j_exc(:,:,xx-yy-zz) + j_exc = 0.0; j_ID = 0 + + if (rank == 0) then + write(6, *) '==> J-values are recieved' + write(6, *) '' + write(6, *) ' ::::::::: J values (meV) ::::::::::::' + write(6, *) '' + write(6, *) ' ion1 <---> ion2 J distance (A)' + end if + do i = 1, no_of_nbd + read(10002, *) nbd_dis(i), similar_bonds + do ii = 1, similar_bonds + + read(10002, *) atom1, atom2, j_value, anisotropy(1:3) + ab = atom1 + call lu(ab(1:1), out, 'U' ) + atom1(1:1) = out + call lu(ab(2:2), out, 'L' ) + atom1(2:2) = out + + ab = atom2 + call lu(ab(1:1), out, 'U' ) + atom2(1:1) = out + call lu(ab(2:2), out, 'L' ) + atom2(2:2) = out + atom = (/ atom1, atom2 /) + + do j = 1, 2 + do k = 1, nspecies + if (atom(j).eq.species(k)) then + j_ID(j) = k + end if + end do + end do + if (rank == 0) write(6, 501) species(j_ID(1)), species(j_ID(2)), j_value, & + nbd_dis(i) +501 format(5X,A2,X,'<--->',X,A2,2X,f8.3,f9.5) + j_exc(i, j_ID(1), j_ID(2), 1:3) = j_value*anisotropy*s(j_ID(1))*s(j_ID(2)) + j_exc(i, j_ID(2), j_ID(1), 1:3) = j_exc(i, j_ID(1), j_ID(2), 1:3) + + end do + end do + + !closing of j_exchange file + close(10002) + + nbd_dis_max = maxval(nbd_dis, dim=1) + + ! Extent of nbd cells around a lattice point + nbd_cell_x = ceiling(nbd_dis_max/lp(1)) + nbd_cell_y = ceiling(nbd_dis_max/lp(2)) + nbd_cell_z = ceiling(nbd_dis_max/lp(3)) + + ! Simulation box start from_ till to_ + fromx = nbd_cell_x + 1; tox = sc(1) + nbd_cell_x + fromy = nbd_cell_y + 1; toy = sc(2) + nbd_cell_y + fromz = nbd_cell_z + 1; toz = sc(3) + nbd_cell_z + + if (rank == 0) then + + write(6, *) '' + write(6, *) ' Note: Kindly check the above listed J values' + write(6, *) ' ````` from j_exchange file before proceeding further!' + write(6, *) '' + + end if + + end subroutine j_values diff --git a/src/lu.f90 b/src/lu.f90 new file mode 100644 index 0000000..cd7d3f0 --- /dev/null +++ b/src/lu.f90 @@ -0,0 +1,42 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 + + !SMALL CAPS <----> LARGE CAPS + subroutine lu(txtR, txtL, case_LU) + + implicit none + + character (len=*), intent(in) :: txtR + character (len(txtR)), intent(out) :: txtL + character(len=54):: s + character, intent(in) :: case_LU + integer :: i, j + + s = ' abcdefghijklmnopqrstuvwxyz ABCDEFGHIJKLMNOPQRSTUVWXYZ' + txtL = txtR + do i = 1, len(txtR) + do j = 1, len(s) + if( txtR(i:i).eq.s(j:j) )then + if( (j.le.27) .and. (case_LU.eq.'U') ) txtL(i:i) = s(j+27:j+27) ! into upper case + if( (j.gt.27) .and. (case_LU.eq.'L') ) txtL(i:i) = s(j-27:j-27) ! into lower case + end if + end do + end do + end subroutine lu + diff --git a/src/mc.f90 b/src/mc.f90 new file mode 100644 index 0000000..f6c2712 --- /dev/null +++ b/src/mc.f90 @@ -0,0 +1,45 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 mc(at_step, acceptance_count_) + + use init + + implicit none + + integer, intent(in) :: at_step + integer, intent(out) :: acceptance_count_ + + ! Perform Monte Carlo + call Monte_Carlo(acceptance_count_) + + after_equilibration_step: if((at_step.gt.tmcs_eq).and.& + (mod(real(at_step), real(to_cal)).eq.0))then + + total_calculations = total_calculations + 1 + eng = 0; net_mag = 0 + call get_tot_energy(eng) + call get_tot_magnetisation(net_mag) + + call evaluate_observables('eng') + call evaluate_observables('mag') + + end if after_equilibration_step + + end subroutine mc diff --git a/src/overrelaxation.f90 b/src/overrelaxation.f90 new file mode 100644 index 0000000..e2ee873 --- /dev/null +++ b/src/overrelaxation.f90 @@ -0,0 +1,77 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 overrelaxation + + use init + + implicit none + + integer :: i, j, k, l, m + + real(8) :: A_ovr(3), Si(3), rn_ovr + + overrelaxation_method : do m = 1, overrelaxed + ! Duplicate Names Between Modules are not Allowed (DNBM) + ! it causes variables to exchange between the modules/threads + call get_random_num(fromx*1d0, tox*1d0, rn_ovr); i = int(rn_ovr) + if(i.lt.fromx) then + i = fromx + elseif(i.gt.tox)then + i = tox + end if + + call get_random_num(fromy*1d0, toy*1d0, rn_ovr); j = int(rn_ovr) !DNBM + if(j.lt.fromy) then + j = fromy + elseif(j.gt.toy)then + j = toy + end if + + call get_random_num(fromz*1d0, toz*1d0, rn_ovr); k = int(rn_ovr) !DNBM + if(k.lt.fromz) then + k = fromz + elseif(k.gt.toz)then + k = toz + end if + + call get_random_num(1d0, lattice_per_unit_cell*1d0, rn_ovr); l = int(rn_ovr) !DNBM + if(l.lt.1) then + l = 1 + elseif(l.gt.(lattice_per_unit_cell))then + l = lattice_per_unit_cell + end if + + A_ovr = 0 + + call get_ovrr_vec(i, j, k, l, A_ovr) + + + + Si(1:3) = ion(1:3, i, j, k, l) + + Si(1:3) = (2* & + (dot_product(A_ovr(1:3), Si(1:3))/dot_product(A_ovr(1:3), A_ovr(1:3)) )* & + A_ovr(1:3)) - Si(1:3) + + ion(1:3, i, j, k, l) = Si(1:3)/sqrt(dot_product(Si(1:3), Si(1:3))) + + end do overrelaxation_method + + end subroutine overrelaxation diff --git a/src/parameters.f90 b/src/parameters.f90 new file mode 100644 index 0000000..cc63b82 --- /dev/null +++ b/src/parameters.f90 @@ -0,0 +1,56 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 parameters + + use init + + implicit none + + if (para) then + lbl = ' KbT/J' + j_exc = j_exc/para_value + kb = 1d0 + mb = 1d0 + g_factor = 1d0 + if(single_ion_anisotropy) sia_vec = sia_vec/para_value + if(Zeeman) h = h/para_value + else + lbl = " K" + ! Converting 'j_exc', 'sia' from meV into eV + j_exc = j_exc*(1d-3) + if(single_ion_anisotropy) sia_vec = sia_vec*(1d-3) + + ! Bhor magneton + mb = 5.7883818060d-5 + if(rank == 0) then + write(6, *) '' + write(6, *) "REMEMBER! For 'J in meV' case." + write(6, *) ' ALL variables should be provided in the terms of milli orders (meV), ' + write(6, *) " in 'input' file. eg., for 1meV put only 1 not 0.001" + write(6, *) '' + write(6, *) ' *******************************' + write(6, *) " KINDLY CHECK THE 'input' FILE" + write(6, *) ' (IGNOR, if it is already DONE!)' + write(6, *) ' *******************************' + write(6, *) '' + end if + end if + + end subroutine parameters diff --git a/src/read_input.f90 b/src/read_input.f90 new file mode 100644 index 0000000..ca7319e --- /dev/null +++ b/src/read_input.f90 @@ -0,0 +1,98 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 read_input + + use init + + implicit none + + integer :: i + character(len=3) :: label + + open(unit=0, file='in.ether', status='old', action='read') + + read(0, *) model ! Ising/Heisenberg model + + ! CONVERTING STRING VARIBLES INTO LOWER CASE + call lu(model, label, 'L' ) !L: lower case, U: upper case + model = label + Ising= .FALSE.; Heisenberg = .FALSE. + if(model.eq.'i') Ising = .TRUE. + if(model.eq.'h') Heisenberg = .TRUE. + + read(0, *) tmcs ! total Monte Carlo Steps per spins + read(0, *) tmcs_eq ! total Monte Carlo Steps per spins for equilibration + read(0, *) ht, lt, tint ! higher temp, lower temp., temp. interval + read(0, *) (s(i), i = 1, nspecies) ! spin-moment + read(0, *) n_speci_incl ! no. of species to include + + allocate(specicies_to_include(n_speci_incl), stgg(n_speci_incl)) + + read(0, *) (specicies_to_include(i), i = 1, n_speci_incl) ! species lable to include in calculation + read(0, *) sc(1), sc(2), sc(3) ! supercell dimension + + if(mod(sc(1), 2).ne.0.or.mod(sc(2), 2).ne.0.or.mod(sc(3), 2).ne.0) then + write(6, *) '' + write(6, *) ' ERROR' + write(6, *) ' ~~~~~' + write(6, *) '' + write(6, *) ' Dimension of supercell along x, y, and z must be in the multiples of 2' + write(6, *) ' STOPPING now' + write(6, *) '' + stop + end if + + read(0, *) staggered ! logic for staggered magnetization + read(0, *) bc(1), bc(2), bc(3) ! boundary condition along x, y, z + + ! CONVERTING STRING VARIBLES INTO LOWER CASE + call lu(bc(1), label, 'L' ) + bc(1) = label + call lu(bc(2), label, 'L' ) + bc(2) = label + call lu(bc(3), label, 'L' ) + bc(3) = label + + read(0, *) sample ! sample per unit MC calculations) + + read(0, *) angle, dphi ! logic, least angle window (helps in convergence) + + dphi = convert_to_rad*dphi + + read(0, *) to_cal ! at this step interval all observables will be calculated + read(0, *) Zeeman, h(1:3) ! Magnetic field logic, Mx, My, Mz (in T) + read(0, *) g_factor ! g-factor + read(0, *) single_ion_anisotropy ! Logic(.T./.F.) + read(0, *) para, para_value ! Logic, parameter value --> calculation in terms of parameter + read(0, *) ovrr, overrelaxed, ovrr_start_interval ! overrelaxed (T/F), overrelaxed step, interval of overrelaxation starts + !read(0, *) EXalgo, exchange_interval, temp_ex !Exchange algo (logic), interval to exchange replicas, for using given equi-spaced temprature (T/F) + !read(0, *) beta_critria, optbeta, beta_file !Beta(M) convergence criteria, optimization steps for beta(M), if Beta.dat file exists(logic) + + if((para_value.eq.0).and.para.and.rank == 0) then + write(6, *) "" + write(6, *) "ERROR:: parameter value is set ON and value & + cannot be zero. Kindly have a look into the input file" + write(6, *) "" + STOP + end if + + close(0) ! closing of input file + + end subroutine read_input diff --git a/src/read_structure.f90 b/src/read_structure.f90 new file mode 100644 index 0000000..0822406 --- /dev/null +++ b/src/read_structure.f90 @@ -0,0 +1,76 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 read_structure(nsp, nions, lp, abc, coordx, coordy, coordz, ions, species_name, rank_) + + implicit none + + integer :: tions, i, j + integer, intent(in) :: nsp, nions, rank_ + integer, intent(out) :: ions(0:nsp) + + real(8), intent(out) :: lp(3), abc(3, 3), coordx(nions), coordy(nions), coordz(nions) + real(8) :: wt + + character(len=200) :: title, label + character(len=2), intent(out) :: species_name(0:nsp) + + logical :: file_present + + inquire(file='structure.vasp', exist=file_present) + if(.not.file_present) then + write(6, *) "==> 'structure.vasp' is not present" + write(6, *) " STOPPING now" + write(6, *) "" + stop + end if + + open(unit=10001, file='structure.vasp', status='old', action='read') + read(10001,*) title + read(10001,*) wt + + lattice_constant : do i = 1,3 + read(10001,*) (abc(i,j), j = 1,3) + enddo lattice_constant + + lp(1) = sqrt(dot_product(abc(1, 1:3), abc(1, 1:3))) ! lattice parameter a + lp(2) = sqrt(dot_product(abc(2, 1:3), abc(2, 1:3))) ! lattice parameter b + lp(3) = sqrt(dot_product(abc(3, 1:3), abc(3, 1:3))) ! lattice parameter c + + read(10001,*) (species_name(i), i = 1,nsp) + species_name(0) = 'X' + ions = 0 + read(10001,*) (ions(i), i =1, nsp) + tions = sum(ions) + + coordx = 0; coordy = 0; coordz = 0 + read(10001, *) label + if(trim(adjustl(label)).eq.'Direct')then + write(6, *) 'Choose POSCAR in cartesian co-ordinate only' + write(6, *) 'STOPPING' + stop + endif + + do i = 1, tions + read(10001,*) coordx(i), coordy(i), coordz(i) !cartesian co-ordinates + end do + close(10001) + if (rank_ == 0) write(6, *) '==> Reading structure is completed!' + + end subroutine read_structure diff --git a/src/startup.f90 b/src/startup.f90 new file mode 100644 index 0000000..c9aef38 --- /dev/null +++ b/src/startup.f90 @@ -0,0 +1,147 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 startup(values) + + use init + + implicit none + + integer,dimension(8), intent(in) :: values + + write(6, *) "" + write(6, *) "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + write(6, *) "~~~~~ ETHER ~~~~~" + write(6, *) "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + write(6, *) "" + write(6, '(" PROGRAM STARTED on date ",i2,"-",i2,"-",i4)') values(3),values(2), values(1) + write(6, "(' at time ',i2,' hrs. ',i2,' min. ',i2,' sec. ')") values(5:7) + write(6, *)"" + + if(Ising) then + write(6, *)"Based on Ising model (H = JSi.Sj)" + elseif(Heisenberg) then + write(6, *)"Based on classical Heisenberg model (H = JSi*Sj)" + else + write(6, *) "WARNING: unable to get right model ID" + write(6, *) "setting model to default 'Heisenberg model (ID = H)'" + model = 'h' + end if + + write(6, *)"" + write(6, *) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + write(6, *) 'Working on Monte Carlo steps' + write(6, *) 'per Spin (MCS)' + write(6, *) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + write(6, "(1X,A21,I8)") 'Total MCS steps: ', tmcs + write(6, "(1X,A21,I8)") 'Equilibration steps: ', tmcs_eq + write(6, *) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + write(6, *) '' + write(6, '(I4,1x,I4,1x,I4," :: Lx, Ly, Lz <==> Supercell size")') sc + + write(6, *)'' + write(6, '(2X,A2,3x,A2,3x,A2,2x"(boundary conditions along x,y, and z-axis, respectively)")') bc(1:3) + write(6, *)" 'o' => open; 'c' => closed" + write(6, *)'' + write(6, *)'==> Spins are set in random configurations' + + if((angle.eqv..FALSE.).and.Heisenberg) then + write(6, *) '==> Spin vectors will be chosen as per George Marsaglia method' + write(6, *)'' + write(6, *) ' > George Marsaglia, The Annals of Mathematical Statistics' + write(6, *) ' Vol. 43, No. 2, 645-646 (1972).' + write(6, *)'' + elseif(angle) then + write(6, '(" ==> +/-",f6.2," deg. angle has been chosen for the convergence")') dphi/convert_to_rad + write(6, *)'' + write(6, *) ' > Application of the Monte Carlo Method in Statistical Physics' + write(6, *) ' by Kurt Binder, Second Edition' + write(6, *)'' + end if + + if(Ising) then + write(6, *) '==> Spin states will be chosen from two states i.e., +S and -S' + end if + + if (ovrr.and..not.Ising)then + write(6, '(" ==> Overrelaxation method with ",I4," OVR steps has been considered along with & + Metropolis algorithm")') overrelaxed + write(6, *)'' + write(6, *) ' > Michael Creutz, Physical Review D Vol. 36, 2 (1987).' + write(6, *) ' > J. L. Alonso and A. Tranacon, Physical Review B Vol. 53, 5 (1996).' + write(6, *) ' > Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenblunth,' + write(6, *) ' Augusta H. Teller, and Edward Teller, The Journal of Chemical Physics' + write(6, *) ' Vol. 21, 1087 (1953).' + write(6, *)'' + else + write(6, *) '==> Metropolis algorithm will be used' + write(6, *)'' + write(6, *) ' > Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenblunth,' + write(6, *) ' Augusta H. Teller, and Edward Teller, The Journal of Chemical Physics' + write(6, *) ' Vol. 21, 1087 (1953).' + write(6, *)'' + end if + + if(EXalgo) then + write(6, *) '==> Exchange alogorithm will be used along with Metropolis algorithm' + write(6, '(" with exchange interval", i4)') exchange_interval + write(6, *) '' + write(6, *) ' > Koji Hukushima and Koji Nemoto, Journal of the Physical Society' + write(6, *) ' of Japan 65, 1604-1608 (1996).' + write(6, *) '' + + if(temp_ex.eqv..FALSE.)then + write(6, '(" ==> Temperature points will be generated as per the method given in the appendix")') + write(6, *) '' + write(6, *) ' > Koji Hukushima, Physical Review E 60, 4 (1999).' + write(6, *)'' + else + write(6, *) '==> Equi-spaced temperature points have been considered for Exchange algorithm' + end if + + end if + + if(staggered)then + write(6, *)'==> Staggered magnetization is SELECTED!' + end if + + if(Zeeman)then + if(para)then + lbl = ' as H/J' + else + lbl = " in Tesla" + end if + write(6, '(" ==> Magnetic field ::", 3f10.5, " (Mx, My, Mz)",A9," is SELECTED!")') h, lbl + end if + + if(single_ion_anisotropy)then + if(para)then + lbl = " as SIA/J" + else + lbl = " in meV" + end if + write(6, '(" ==> SIA (Single Ion Anisotropy) ", A9," is SELECTED!")') lbl + end if + + if(para)then + write(6, *) '==> Calculations in terms of parameters is ON!' + write(6, '(" ==> Parameters w.r.t. ",f8.3, " meV has been taken")') para_value + end if + + end subroutine startup diff --git a/src/update_spin_details.f90 b/src/update_spin_details.f90 new file mode 100644 index 0000000..d7f59b7 --- /dev/null +++ b/src/update_spin_details.f90 @@ -0,0 +1,53 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 update_spin_details(ref1, ref2) + + use init + + implicit none + + real(8) :: rnum, sx, sy, sz, phi + real(8), intent(in) :: ref1(5) + real(8), intent(out) :: ref2(5) + + ref2 = ref1 ! storing species no. + if(Ising) then + ref2(3) = -ref1(3) + else + if(angle.eqv..FALSE.) then + call George_Marsaglia(sx, sy, sz) + ref2(1) = sx; ref2(2) = sy; ref2(3) = sz + ref2(5) = 0 + else + ref2(5) = ref1(5) ! old phi + call get_random_num(-1d0, +1d0, rnum) + ref2(3) = rnum ! cos(theta) : S(z) + call get_random_num(-1d0, +1d0, rnum) + phi = ref2(5) + (dphi*rnum) ! older_phi + delta*rnum + ref2(1) = sqrt(1d0-(ref2(3)**2))*cos(phi) ! S(x) + ref2(2) = sqrt(1d0-(ref2(3)**2))*sin(phi) ! S(y) + ref2(5) = phi ! storing phi value + end if + end if + + ! Normalising to mod 1 + ref2(1:3) = ref2(1:3)/sqrt(dot_product(ref2(1:3), ref2(1:3))) + + end subroutine update_spin_details diff --git a/src/write_gss.f90 b/src/write_gss.f90 new file mode 100644 index 0000000..ba41cb0 --- /dev/null +++ b/src/write_gss.f90 @@ -0,0 +1,89 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 write_gss(tempi) + + use init + + implicit none + + integer :: i, j, k, l, m + integer, intent(in) :: tempi + gss_ID = 10009 + + if(tempi.le.0) then + + open(unit=gss_ID, file='gss.dat', status='replace', action='write') + + if(para)then + lbl = ' KbT/J' + else + lbl = " K" + end if + + ! STORING crystal information in GSS (Ground Spin States) + write(gss_ID, '(I5,1x, I5, 1x, I5, 1x, I5)') nscan, sc + write(gss_ID, '(I5,1x, I5, 1x, I5)') nbd_cell_x, nbd_cell_y, nbd_cell_z + write(gss_ID, '(I5,1x, I5, 1x, I5, 1x, I5, 1x, I5, 1x, I5)') fromx, & + fromy, fromz, tox, toy, toz + write(gss_ID, '(I5,1x, I5, f11.5, 1x, f11.5, 1x, f11.5,1x, A6)') & + lattice_per_unit_cell, nspecies, ht, lt, tint + write(gss_ID, '(A20)') lbl + write(gss_ID, '(10I5)') ionn(1:) + do i = 1,3 + write(gss_ID,'(3f10.5)') (abc(i,j), j= 1,3) + enddo + write(gss_ID,*) '# lattice points' + do i = 1, sc(1) + 2*nbd_cell_x + do j = 1, sc(2) + 2*nbd_cell_y + do k = 1, sc(3) + 2*nbd_cell_z + do l = 1, lattice_per_unit_cell + + write(gss_ID,'(f11.6, 1x, f11.6, 1x, f11.6, 1X, i8)') & + ion(6:8, i, j, k, l), int(ion(0, i, j, k, l)) + + end do + end do + end do + end do + return + + end if + + !GSS.dat + m = 0 + write(gss_ID, '(f11.5)') temp_T(tempi) + do i = fromx, tox + do j = fromy, toy + do k = fromz, toz + do l = 1, lattice_per_unit_cell + + write(gss_ID, '(f11.6, 1x, f11.6, 1x, f11.6,1X, i8, 1x, i8)') & + global_spn(m + 1 + lspn), global_spn(m + 2 + lspn), & + global_spn(m + 3 + lspn), int(global_spn(m + 4 + lspn)) + m = m + 4 + + end do + end do + end do + end do + + if(tempi.eq.nscan) close(gss_ID) + + end subroutine write_gss diff --git a/src/write_initial_conf.f90 b/src/write_initial_conf.f90 new file mode 100644 index 0000000..8063dcb --- /dev/null +++ b/src/write_initial_conf.f90 @@ -0,0 +1,50 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 + + ! WRITING INITIAL STRUCTURE FILE + subroutine write_initial_conf + + use init + + implicit none + + integer :: i, j, k, l + + open(unit=10003, file='starting_spin_conf.xsf', status='unknown') + write(10003, *) 'ATOM' + do i = 1, sc(2) + 2*nbd_cell_x + do j = 1, sc(2) + 2*nbd_cell_y + do k = 1, sc(3) + 2*nbd_cell_z + do l = 1, lattice_per_unit_cell + + write(10003,10021) species(int(ion(4, i, j, k, l))), & + ion(6:8, i, j, k, l), ion(1:3, i, j, k, l) + + end do + end do + end do + end do + close(10003) + +10031 format(12A7) +10021 format(A5,7f13.7) + + write(6, *) "==> Initial spin states have be written into 'starting_spin_conf.dat'" + + end subroutine write_initial_conf diff --git a/src/write_nbd.f90 b/src/write_nbd.f90 new file mode 100644 index 0000000..2b0b4a3 --- /dev/null +++ b/src/write_nbd.f90 @@ -0,0 +1,64 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 write_nbd + + use init + + implicit none + + integer :: i, j, k, l, m + + open(unit=10006, file='nbd.dat', status='unknown') + write(6, *) "==> Neighbourhood informations are writing into the 'nbd.dat' file" + + write(10006, *) "This file can be used to confirm the neighbourhood details of any central ion." + write(10006, *) "To do so, use 'XCrySDen' (http://www.xcrysden.org/) and open the generated" + write(10006, *) "'starting_spin_conf.xsf' as file > Open Structure > Open XSF > starting_spin_conf.dat." + write(10006, *) "Check nbd's details by clicking the (at bottom) and see the 'Selected Atom No.' informations. " + write(10006, *) '' + + nbd_inf : do i = 1, sc(1) + 2*nbd_cell_x + do j = 1, sc(2) + 2*nbd_cell_y + do k = 1, sc(3) + 2*nbd_cell_z + do l = 1, lattice_per_unit_cell + + write(10006,10) int(ion(0, i, j, k, l)) + write(10006,*)"~~~~~~~~~~~~~~" + + nbd_write : do m = 1, no_of_nbd + write(10006, 11) m, & + nn(m, int(ion(0, i, j, k, l)), & + 1:nn(m, int(ion(0, i, j, k, l)), 0, 0), & + 0) + +10 format(" ION no. ",i5) +11 format(" For bond length no. ",i2," nbds are ==> "20i7) + end do nbd_write + + write(10006,*)'' + + end do + end do + end do + end do nbd_inf + + close(10006) + + end subroutine write_nbd diff --git a/src/write_output_files.f90 b/src/write_output_files.f90 new file mode 100644 index 0000000..90587cd --- /dev/null +++ b/src/write_output_files.f90 @@ -0,0 +1,83 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 write_output_files(tempi) + + use init + + implicit none + + integer :: i + integer, intent(in) :: tempi + + if(tempi.lt.1) then + + open(unit=20001, file='magnetization.dat', status='replace',action='write') + open(unit=20002, file='energy.dat', status='replace', action='write') + open(unit=20003, file='moment.dat', status='replace', action='write') + + m_head(1) ='# 1 Temp.' + m_head(2) ='2 |M|'; m_head(3)='3 χ' + m_head(4) ='4 Δ|M|'; m_head(5) ='5 Δχ' + m_head(6) ='6 U(|M|)'; m_head(7) ='7 ΔU(|M|)' + + e_head(1) ='# 1 Temp.' + e_head(2) ='2 E' + + if(para)then + e_head(3)='3 Cv' ; e_head(5) ='5 ΔCv' + else + e_head(3)='3 Cv/kB' ; e_head(5) ='5 Δ(Cv/kB)' + end if + + e_head(4) ='4 ΔE' + e_head(6) ='6 U(E)'; e_head(7) ='7 ΔU(E)' + e_head(8) ='8 Accpt. (%)' + + write(20001,103) (m_head(i), i = 1, 7) + write(20002,104) (e_head(i), i = 1, 8) +103 format(A13,A13,A13,1x,A13,1x,A13,2x,A13,A13) +104 format(A13,A13,A13,A13,1x,A13,1x,A13,A13,A13) + + write(20003, *) '# spN = Nth species, M = magnitude of magnetic moment vector (Mx, My, Mz)' + write(20003, *) '# Temp, sp1(Mx,My,Mz), sp1(Mx,My,Mz),..spN(Mx,My,Mz), sp1|M|, sp2|M|,..spN|M|' + + return + end if + + !MAGNETIC MOMENT + write(20003, 101) temp_T(tempi), (mm_vector_avg_T(i, 1:3, tempi)& + /(product(sc)*ionn(i)), i = 1, nspecies), & + (sqrt(dot_product(mm_vector_avg_T(i, 1:3, tempi), mm_vector_avg_T(i, 1:3, tempi)))& + /(product(sc)*ionn(i)), i = 1, nspecies) + + !MAGNETIC + write(20001,102) temp_T(tempi), s_mag_avg_T(tempi), & + s_chi_T(tempi), err_mag_avg_T(tempi), & + err_chi_T(tempi), s_U_mag_T(tempi), err_U_mag_T(tempi) + + !ENERGY + write(20002,102) temp_T(tempi), s_eng_avg_T(tempi), & + s_cv_T(tempi), err_eng_avg_T(tempi), & + err_cv_T(tempi), s_U_eng_T(tempi), err_U_eng_T(tempi), acceptance_ratio(tempi) + +102 format(f11.5, 1x, es12.5, 1x,es12.5, 1x,es12.5, 1x,es12.5, 1x,es12.5, 1x,es12.5, 1x, f7.2) +101 format(f11.5, 2X,500f9.3) + + end subroutine write_output_files diff --git a/src/write_spins_at_K.f90 b/src/write_spins_at_K.f90 new file mode 100644 index 0000000..b3f6997 --- /dev/null +++ b/src/write_spins_at_K.f90 @@ -0,0 +1,77 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 write_spins_at_K(tempi) + + use init + + implicit none + + integer, intent(in) :: tempi + integer :: i, j, k, l, sce(3) + + if(initiate_spin_files) then + + initiate_spin_files = .FALSE. + !FORMATION OF FILES @ temp. K + spin_file_ID = 2001 + tempi + write(filename, 100) temp +100 format( f9.4,'spK.xsf') + open(file = adjustl(filename), unit = spin_file_ID) + write(spin_file_ID, *) 'CRYSTAL' + write(spin_file_ID, *) 'PRIMVEC' + sce = sc ! Supercell extent (sce) + if(bc(1).eq.'o') sce(1) = sc(1)+1 + if(bc(2).eq.'o') sce(2) = sc(2)+1 + if(bc(3).eq.'o') sce(3) = sc(3)+1 + write(spin_file_ID, *) (abc(1,j)*(sce(1)), j= 1,3) + write(spin_file_ID, *) (abc(2,j)*(sce(2)), j= 1,3) + write(spin_file_ID, *) (abc(3,j)*(sce(3)), j= 1,3) + write(spin_file_ID, *) 'CONVEC' + write(spin_file_ID, *) (abc(1,j)*(sce(1)), j= 1,3) + write(spin_file_ID, *) (abc(2,j)*(sce(2)), j= 1,3) + write(spin_file_ID, *) (abc(3,j)*(sce(3)), j= 1,3) + write(spin_file_ID, *) 'PRIMCOORD' + write(spin_file_ID, *) product(sc)*lattice_per_unit_cell, " 1" + return + + end if + + ! SPINS + do i = fromx, tox + do j = fromy, toy + do k = fromz, toz + do l = 1, lattice_per_unit_cell + + if (ion(4, i, j, k, l).ne.0) then + + write(spin_file_ID, 101) species(int(ion(4, i, j, k, l))), & + ion(6:8, i, j, k, l), ion(1:3, i, j, k, l) + + end if + + end do + end do + end do + end do +101 format(A5,7f13.7) + + close(spin_file_ID) + + end subroutine write_spins_at_K diff --git a/src/zeroes.f90 b/src/zeroes.f90 new file mode 100644 index 0000000..5bc4972 --- /dev/null +++ b/src/zeroes.f90 @@ -0,0 +1,61 @@ +! 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 (msharma1@ph.iitr.ac.in) + +! 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 zeroes(zero_opt) + + use init + + implicit none + + character(len=*), intent(in) :: zero_opt + + select case(zero_opt) + + case('sample') + + s_eng_avg = 0.0; e_eng_avg = 0.0 + s_eng2_avg = 0.0; e_eng2_avg = 0.0 + s_U_eng = 0.0; e_U_eng = 0.0 + s_cv = 0.0; e_cv = 0.0 + + s_mag_avg = 0.0; e_mag_avg = 0.0 + s_mag2_avg = 0.0; e_mag2_avg = 0.0 + s_U_mag = 0.0; e_U_mag = 0.0 + s_chi = 0.0; e_chi = 0.0 + + case('eng_mag') + + eng_avg = 0.0; eng2_avg = 0.0; eng4_avg = 0.0 + mag_avg = 0.0; mag2_avg = 0.0; mag4_avg = 0.0 + total_calculations = 0 + + case default + + if (rank == 0) then + write(6, *) '' + write(6, "(' ==> Found unknown case tag :',A8 )") zero_opt + write(6, *) "in 'zeroes' subroutine" + write(6, *) 'STOPPING now' + stop + end if + + end select + + + end subroutine zeroes