784 lines
26 KiB
Fortran
Executable File
784 lines
26 KiB
Fortran
Executable File
! dynamics.f90
|
|
!
|
|
!
|
|
!****************************************************************************
|
|
!
|
|
! PROGRAM: dynamics
|
|
!
|
|
! PURPOSE: Simulation of the asymmetry of an artificial spinglass.
|
|
! DYNAMICS assumes the spinglass to have an FCC lattice.
|
|
! The dimensions of the lattice are w * w * d, along the
|
|
! x-, y- and z-axis respectively. The magnetic moments
|
|
! are randomly distributed over the latticepoints, the muons
|
|
! are placed on the centers of the FCC-cube.
|
|
! The directions of the magnetic moments is choosen randomly
|
|
! over the whole sphere.
|
|
! The program calculates the magnetic field at the site of
|
|
! muon by adding all dilopar contributions from about 300
|
|
! magnetic moments which nearest by the muonsite. Periodic
|
|
! boundary conditions are applied in the x- and y-direction.
|
|
! The z-direction is assumed to perpendicular to a thin
|
|
! film surface.
|
|
! The dynamics of the magnetic spins is included in one of
|
|
! the following ways:
|
|
! For fluctuationrates larger then 100 MHz,
|
|
! a timestep tau is choosen from a
|
|
! log distrubution ( tau = - ln(random) / fluctuationrate )
|
|
! The muon then rotates for tau microseconds, after all spins
|
|
! are rotated over an angle between - dtetha en dtheta and - dphi and dphi.
|
|
! This process is repeated until the total time is 10 microsecods or more.
|
|
! Output of the muon position is done about every time_resolution microsecond.
|
|
! For fluctuationrates smaller then 100 MHz
|
|
! the muons rotate 1000 times for time_resolution microsecond, after each rotation
|
|
! a fraction (= fluctuationrate / 100) of the magnetic ions are rotated
|
|
! over an angle between - dtetha en dtheta and - dphi and dphi.
|
|
! After each fluctuation the fields at the muonsites are recalculated.
|
|
! "deporization" functions are calculated for
|
|
! left-right, up-down and forward-backward detectors,
|
|
! being the x-, y- and z-components of the muon spin vector.
|
|
! For arbitrary direction one has to take the scalar product of
|
|
! that specific direction with the results produced by this program
|
|
!
|
|
! USE: The parameters used for the simulation are supposed to be on
|
|
! file with the generic name <calculation>.inp.
|
|
! The program can be started in two ways:
|
|
!
|
|
! typing DYNAMICS
|
|
! the user will be prompted for the name of the calculation
|
|
!
|
|
! typing DYNAMICS <calculation>
|
|
! the name of the calculation will be read from the commandline.
|
|
!
|
|
! Output will be written on <calculation>.out and on separate files
|
|
! (for each set of parameters) named <calculation>_###.g_t, where
|
|
! ### can a unique number according to the following rules:
|
|
! If a file \simulations\counter.his can be opened, the program will
|
|
! the number in this file and uses that as a start for numbering
|
|
! the *.g_t files. The program will update \simulations\counter.his.
|
|
! If that file is not present, the program will start at number 1.
|
|
!
|
|
! INPUT: For each simulation the following set of parameters has to be
|
|
! given on one line in the file <calculation>.inp
|
|
!
|
|
! lattice parameter [nm]
|
|
! magnetic moment [Bohr-magneton]
|
|
! external field, three component [tesla]
|
|
! thickness d [nm]
|
|
! width w [nm]
|
|
! concentration [at.%]
|
|
! number of muons #
|
|
! initial muon spin direction in
|
|
! spherical coordinates, theta, phi [degree]
|
|
! note that the z-axis is perpendicular to the film
|
|
! muon stopping range, from d1 to d2 [nm]
|
|
! d1 and d2 are note restricted by 0 and d, stopping outside
|
|
! the actual sample is possible
|
|
! fluctuationrate [ 1/ microsecond ]
|
|
! fluctuation amplitude, in
|
|
! spherical coordinates, d-theta, d-phi [degree]
|
|
!
|
|
! Lines with parameters can be interlaced with comments,
|
|
! Commentlines should have a ! at position 1.
|
|
!
|
|
!
|
|
!
|
|
!****************************************************************************
|
|
|
|
program dynamics
|
|
|
|
Use DFPORT
|
|
Use DFLIB
|
|
|
|
implicit none
|
|
|
|
! Variables
|
|
|
|
integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
|
|
& max_muons = 10000, & ! maximum number of muons
|
|
& max_nn = 500, & ! maximum number of nearest neighbours
|
|
& n_time_steps = 500 ! number of time steps calculated
|
|
! Should be future variable
|
|
! Should be future variable
|
|
|
|
|
|
! Structure to store the position (as lattice site-indexes)
|
|
! and the direction-cosines of each spin.
|
|
|
|
structure /spin/
|
|
integer*4 x,y,z
|
|
real*8 theta,phi,dir(3)
|
|
end structure
|
|
|
|
structure /muon/
|
|
integer*4 x,y,z, ns, s(max_nn)
|
|
real*8 dir(3), r(3,max_nn), r_2(max_nn), r_5(max_nn), omega(3)
|
|
end structure
|
|
|
|
! Declarations, maximumnumber of spins: max_spins, maxd is the maximum number of
|
|
! unitcell-distance for which the spin in included in the calculation
|
|
|
|
real*8, time_resolution=0.01 ! time resolution
|
|
! (approximate time between points
|
|
! on the *.g_t file)
|
|
|
|
|
|
character*10 dddd, tttt, zone
|
|
character*4 file_index
|
|
integer*4 dt(8), ifile, l_calc, n_steps, i_step
|
|
character*80 calculation
|
|
character*127 line
|
|
logical unique
|
|
integer*4 i,j,k,l,nw,nd,nsp,n_spin, n_site
|
|
integer*4 iseed, nd1, nd2
|
|
record /spin/ s(max_spins)
|
|
record /muon/ m(max_muons)
|
|
real*8 dummy, a, d, concentration, w, depth1, depth2
|
|
real*8 factor, moment, b_ext(3)
|
|
real*8 fluctuationrate, tau, dphi, dtheta, fraction
|
|
real*8 g_t(3), omega, b_abs, b_sq, ca_sq, his, radiussq
|
|
real*8 t_ini, p_ini, emu(3), Pi
|
|
real*8 step, exp_time, write_time, f_c
|
|
|
|
! Body of dynamics
|
|
|
|
! Read the parameters from the input file
|
|
! with name : <calculation>.inp
|
|
! The output will go to <calculation>.out
|
|
! and <calculation>_###.g_t
|
|
|
|
Write(6,*) ' '
|
|
Write(6,*) ' ---------------------------------------------------------------------'
|
|
Write(6,*) ' | Program field-calculation of muons due to random dynamic spins |'
|
|
Write(6,*) ' | |'
|
|
Write(6,*) ' | Version of November 16 2005 |'
|
|
Write(6,*) ' | |'
|
|
Write(6,*) ' | Input is read from an input file that should be named |'
|
|
Write(6,*) ' | <calculation>.inp and contains for each simulation on |'
|
|
Write(6,*) ' | one line: |'
|
|
Write(6,*) ' | |'
|
|
Write(6,*) ' | lattice-constant [nm], magnetic moment [mu_B], |'
|
|
Write(6,*) ' | ext. field(3) ,thickness, width, c, number_of_muons, |'
|
|
Write(6,*) ' | initial-muon-direction(theta, phi)[degrees], |'
|
|
Write(6,*) ' | (muon-positions from) depth1, (to) depth2 [nm], |'
|
|
Write(6,*) ' | fluctions rate [inverse microsec], |'
|
|
Write(6,*) ' | fluctuation amplitude parallel to film [0..360degr.],|'
|
|
Write(6,*) ' | fluctuation amplitude perpendicular to film |'
|
|
Write(6,*) ' | [0..180degr.] |'
|
|
Write(6,*) ' | |'
|
|
Write(6,*) ' | Lines with a ! in the first position are treated as comments. |'
|
|
Write(6,*) ' | |'
|
|
Write(6,*) ' | <calculation> can be issued as a commandline parameter |'
|
|
Write(6,*) ' ---------------------------------------------------------------------'
|
|
|
|
! files :
|
|
|
|
inquire( file='\simulations\counter.his', exist = unique )
|
|
IF ( unique ) THEN
|
|
open(9,file='\simulations\counter.his',status='old',err=994)
|
|
read(9,*) ifile ! initialize outputfile counter
|
|
ELSE
|
|
ifile = 0
|
|
END IF
|
|
|
|
IF ( iargc() .GT. 0 ) THEN
|
|
call getarg(1, calculation)
|
|
Write(6,*) ' Calculation taken from commandline > ',calculation
|
|
ELSE
|
|
200 write(6,201)
|
|
201 format(' '/' Give name of the calculation > ', \)
|
|
read(5,'(a60)') calculation
|
|
END IF
|
|
|
|
l_calc = index( calculation, ' ') - 1
|
|
|
|
IF ( l_calc .GT. 0 ) THEN
|
|
open(1,file=calculation(1:l_calc)//'.inp',status='old',action='read',err=995 )
|
|
open(2,file=calculation(1:l_calc)//'.out',status='unknown',action='write',err=996)
|
|
|
|
|
|
END IF
|
|
|
|
! initialization of randomumber generator and Pi
|
|
|
|
iseed = 1234567
|
|
Pi = acos( -1.0D+00 )
|
|
|
|
! Read everything from the input file, one line per calculation
|
|
|
|
DO WHILE ( .NOT. Eof(1) ) ! WHILE LOOP(1) over the input file
|
|
|
|
10 read(1,'(a127)',end=999) line
|
|
|
|
IF ( line(1:1) .EQ. '!' ) THEN
|
|
Write(2,'(a)') line
|
|
GOTO 10
|
|
END IF
|
|
|
|
read(line,*,err=998,end=999) a, moment, b_ext, d, w, concentration, &
|
|
& n_site, t_ini, p_ini, depth1, depth2, &
|
|
& fluctuationrate, dphi, dtheta
|
|
|
|
IF ( n_site .GT. 0.8 * max_muons ) n_site = 0.8 * max_muons
|
|
|
|
! Estimate optimum time_step
|
|
|
|
b_abs = sqrt( sum( b_ext * b_ext ) )
|
|
f_c = 135.5 * b_abs + 0.3 * moment * comcentration / (a*a*a)
|
|
|
|
time_step = min( 0.01, 0.1 / f_c )
|
|
|
|
! Initialize randomnumber generator "randomly"
|
|
|
|
call date_and_time( dddd, tttt, zone, dt )
|
|
DO i = 1, dt(8) ! number milliseconds on the clock
|
|
dummy = ran(iseed)
|
|
END DO
|
|
!
|
|
|
|
write(2,100) calculation(1:73),(dt(j),j=1,3),(dt(j),j=5,8)
|
|
100 format(' '/' ',73('-')/' ',a73/' ',73('-')/ &
|
|
& ' Calculation started ',i5,'-',i2,'-',i2, &
|
|
& ' at ',2(i2,':'),i2,'.',i3/' ',73('-')/' ')
|
|
write(2,'(a,f8.3)') ' lattice parameter = ', a
|
|
write(2,'(a,f8.3)') ' magnetic moment = ', moment
|
|
write(2,'(a,3f8.3)') ' external field = ', b_ext
|
|
write(2,'(a,f8.3)') ' concentration = ', concentration
|
|
write(2,'(a,2f8.3)') ' init.muon theta,phi = ', t_ini, p_ini
|
|
write(2,'(a,f8.3)') ' fluctuationrate = ', fluctuationrate
|
|
write(2,'(a,2f8.3)') ' fluctuation amp. = ', dphi, dtheta
|
|
|
|
emu(1) = sin(t_ini*Pi/180.0) * cos(p_ini*Pi/180.0)
|
|
emu(2) = sin(t_ini*Pi/180.0) * sin(p_ini*Pi/180.0)
|
|
emu(3) = cos(t_ini*Pi/180.0)
|
|
|
|
DO j = 1, max_muons
|
|
m(j).dir = emu
|
|
END DO
|
|
|
|
exp_time = 0.0D+00
|
|
write_time = 0.0D+00
|
|
|
|
! update file index for
|
|
|
|
ifile = ifile + 1 ! increase outputfile number
|
|
|
|
IF ( unique ) THEN
|
|
rewind(9)
|
|
write(9,*) ifile ! store for next program
|
|
END IF
|
|
|
|
write(file_index,'(''_'',i3)') ifile ! generate file_name
|
|
|
|
DO j = 2, 4
|
|
IF (file_index(j:j) .EQ. ' ' ) file_index(j:j) = '0'
|
|
END DO
|
|
|
|
open(3,file=calculation(1:l_calc)//file_index//'.g_t', &
|
|
& status='unknown',action='write', err=997 )
|
|
|
|
! output time=0 asymmetries
|
|
|
|
write(3,'(4F19.6)' ) exp_time, emu
|
|
write_time = write_time + time_resolution
|
|
|
|
! Initialize randomnumber generator "randomly"
|
|
|
|
call date_and_time( dddd, tttt, zone, dt )
|
|
|
|
DO i = 1, dt(8) ! number milliseconds on the clock
|
|
dummy = ran(iseed)
|
|
END DO
|
|
|
|
|
|
! make lattice, spinglass, choose muon sites and calculates "interaction matrix"
|
|
|
|
CALL lattice( iseed, d, w, a, concentration, n_spin, nd, nw, s, &
|
|
& n_site, depth1, depth2, m )
|
|
|
|
|
|
write(2,'(a,f8.3)') ' thickness = ', d
|
|
write(2,'(a,f8.3)') ' width = ', w
|
|
write(2,'(a,i10)') ' number of spins = ', n_spin
|
|
write(2,'(a,i10)') ' number of muons = ', n_site
|
|
write(2,'(a,2f8.3)') ' muons penetrate betw. ', depth1, depth2
|
|
write(2,'(a,a)') ' Output will be written on ', &
|
|
& calculation(1:l_calc)//file_index//'.g_t'
|
|
|
|
write(6,'(a,f8.3)') ' lattice parameter = ', a
|
|
write(6,'(a,f8.3)') ' magnetic moment = ', moment
|
|
write(6,'(a,3f8.3)') ' external field = ', b_ext
|
|
write(6,'(a,f8.3)') ' concentration = ', concentration
|
|
write(6,'(a,2f8.3)') ' init.muon theta,phi = ', t_ini, p_ini
|
|
write(6,'(a,f8.3)') ' fluctuationrate = ', fluctuationrate
|
|
write(6,'(a,2f8.3)') ' fluctuation amp. = ', dphi, dtheta
|
|
write(6,'(a,f8.3)') ' thickness = ', d
|
|
write(6,'(a,f8.3)') ' width = ', w
|
|
write(6,'(a,i10)') ' number of spins = ', n_spin
|
|
write(6,'(a,i10)') ' number of muons = ', n_site
|
|
write(6,'(a,2f8.3)') ' muons penetrate betw. ', depth1, depth2
|
|
write(6,'(a,a)') ' Output will be written on ', &
|
|
& calculation(1:l_calc)//file_index//'.g_t'
|
|
|
|
|
|
|
|
! The fluctuations are incorporated as follows:
|
|
|
|
! for rates larger then 1/time_resolution MHz,
|
|
! a timestep tau is choosen from a
|
|
! log distrubution ( tau = - ln(random) / fluctuationrate )
|
|
! The muon then rotates for tau microseconds, after all spins
|
|
! are rotated over an angle between - dtetha en dtheta and - dphi and dphi.
|
|
! This process is repeated until the total time is n_time_steps*time_resolution
|
|
! microsecods or more.
|
|
! Output of the muon position is about every time_resolution microsecond.
|
|
|
|
! for rates smaller then 1/time_resolution MHz
|
|
! the muons rotate 1000 times for time_resolution microsecond, after each rotation
|
|
! a fraction (= fluctuationrate *time_resolution) of the magnetic ions are rotated
|
|
! over an angle between - dtetha en dtheta and - dphi and dphi.
|
|
|
|
! After each fluctuation the fields at the muonsites are recalculated.
|
|
|
|
! "deporization" functions are calculated for
|
|
! left-right, up-down and forward-backward detectors,
|
|
! being the x-, y- and z-components of the muon spin vector.
|
|
|
|
! For arbitrary direction one has to take the scalar product of
|
|
! that specific direction with the results produced by this program
|
|
|
|
|
|
IF ( fluctuationrate .GT. 1.0/time_resolution ) THEN ! Rapid fluctuations
|
|
|
|
fraction = 1.0
|
|
! Start of WHILE loop(2) over exp_time
|
|
|
|
DO WHILE ( exp_time .LT. time_resolution * float(n_time_steps) )
|
|
|
|
tau = - log( ran(iseed) ) / fluctuationrate
|
|
|
|
IF ( exp_time + tau .GT. time_resolution * float(n_time_steps) ) &
|
|
& tau = time_resolution * float(n_time_steps) - exp_time
|
|
|
|
! take at least time_resolution microsec. steps, even if tau is larger
|
|
|
|
n_steps = floor( tau / time_resolution ) + 1
|
|
step = tau / float( n_steps )
|
|
|
|
CALL fields( a, moment, b_ext, s, n_site, m)
|
|
|
|
DO i_step = 1, n_steps
|
|
exp_time = exp_time + step
|
|
call muonrotation( n_site, m, step )
|
|
g_t = 0.0
|
|
DO j = 1, n_site
|
|
DO k = 1, 3
|
|
g_t(k) = g_t(k) + m(j).dir(k)
|
|
END DO
|
|
END DO
|
|
|
|
g_t = g_t / float(n_site)
|
|
|
|
IF ( exp_time .GT. write_time ) THEN
|
|
write(3,'(4F19.6)' ) exp_time, g_t
|
|
write_time = exp_time + time_resolution
|
|
END IF
|
|
|
|
END DO
|
|
|
|
! after tau, change spin directions and repeat the above.
|
|
! however, stop when 10 microsec has been reached.
|
|
|
|
CALL fluctuation( iseed, n_spin, s, dtheta, dphi, fraction )
|
|
|
|
END DO ! END of WHILE loop(2)
|
|
|
|
ELSE ! fluctuationrate < 1/time_resolution
|
|
|
|
fraction = fluctuationrate * time_resolution
|
|
step = time_resolution
|
|
n_steps = n_time_steps
|
|
|
|
DO i_step = 1, n_steps
|
|
exp_time = exp_time + step
|
|
CALL fields( a, moment, b_ext, s, n_site, m)
|
|
CALL muonrotation( n_site, m, step )
|
|
|
|
g_t = 0.0
|
|
|
|
DO j = 1, n_site
|
|
DO k = 1, 3
|
|
g_t(k) = g_t(k) + m(j).dir(k)
|
|
END DO
|
|
END DO
|
|
|
|
g_t = g_t / float(n_site)
|
|
|
|
write(3,'(4F19.6)' ) exp_time, g_t
|
|
CALL fluctuation( iseed, n_spin, s, dtheta, dphi, fraction )
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END DO ! END of WHILE loop(1)
|
|
|
|
STOP ' Program DYNAMICS stopped where it should stop '
|
|
|
|
994 STOP ' FATAL: Cannot open counter.his '
|
|
|
|
995 STOP ' FATAL: Cannot open input file '
|
|
|
|
996 STOP ' FATAL: Cannot open output file '
|
|
|
|
997 Write(2,*) ' Cannot open g_t file '
|
|
STOP ' FATAL: Cannot open g_t file '
|
|
|
|
998 Write(2,*) ' Error in input file '
|
|
STOP ' FATAL: Due to error in input file '
|
|
|
|
999 Write(2,*) ' End of input-file '
|
|
STOP ' STOP End of input file '
|
|
|
|
end program dynamics
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
! LATTICE calculates the actual dimensions of the sample, places magnetic spins
|
|
! randomly according to concentration, gives the spins a random direction
|
|
! in space. It also generates a table of muonsites.
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
Subroutine lattice( iseed, d, w, a, concentration, n_spin, nd, nw, s, &
|
|
& n_site, depth1, depth2, m )
|
|
|
|
! Structure to store the position (as lattice site-indexes)
|
|
! and the direction-cosines of each spin and muon.
|
|
|
|
implicit none
|
|
|
|
integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
|
|
& max_muons = 10000, & ! maximum number of muons
|
|
& max_nn = 500 ! maximum number of nearest neighbours
|
|
|
|
|
|
structure /spin/
|
|
integer*4 x,y,z
|
|
real*8 theta,phi,dir(3)
|
|
end structure
|
|
|
|
structure /muon/
|
|
integer*4 x,y,z, ns, s(max_nn)
|
|
real*8 dir(3), r(3,max_nn), r_2(max_nn), r_5(max_nn), omega(3)
|
|
end structure
|
|
|
|
real*8 d, w, a, concentration, c, depth1, depth2, fraction, radiussq
|
|
real*8 Pi, r(3), r_2, r_3, r_5, help
|
|
integer*4 iseed, nd, nw, nat, n_spin, n_site, nd1, nd2
|
|
integer*4 i, j, k, l, hw, kw, ns
|
|
record /spin/ s(*)
|
|
record /muon/ m(*)
|
|
|
|
Pi = acos( -1.0D+00 )
|
|
c = concentration / 100.0
|
|
|
|
! Calculate the 'rounded' number of spins for a lattice m*m*n for
|
|
! the given concentration.
|
|
! n is the number of atoms (half unitcells) perpendicular
|
|
! to the layer (== z-direction).
|
|
! m is the size of the layer in the x- and y-direction
|
|
|
|
nd = floor(2.0 * d / a ) + 2
|
|
nw = floor(2.0 * w / a ) + 2
|
|
nat = nd * nw * nw / 2
|
|
n_spin = floor( nat * c )
|
|
d = float(nd) * a / 2.0
|
|
w = float(nw) * a / 2.0
|
|
hw = nw / 2
|
|
|
|
IF ( c .GT. 0.0 ) THEN
|
|
radiussq = (( 1.6 * float(max_nn) / c ) / ( 4.0 * Pi / 3.0 ))**(2.0/3.0)
|
|
ELSE
|
|
radiussq = 1D+10
|
|
END IF
|
|
|
|
write(2,*) ' radius = ', sqrt( radiussq )
|
|
|
|
nd1 = floor( 2.0 * depth1 / a )
|
|
nd2 = floor( 2.0 * depth2 / a )
|
|
IF ( mod( nd1 , 2 ) .EQ. 0 ) nd1 = nd1 + 1 ! nd1 should be odd
|
|
IF ( nd2 .LT. nd1 + 1 ) nd2 = nd1 + 1
|
|
depth1 = float(nd1) * a / 2.0
|
|
depth2 = float(nd2) * a / 2.0
|
|
fraction = float(n_site) / (float((nd2-nd1)*nw*nw) / 8.0)
|
|
|
|
|
|
! Place the spins randomly on the fcc-lattice
|
|
! Run over a whole simple cubic lattice in steps
|
|
! of half of the fcc-unitcell.
|
|
! Then take care of the fcc-structure and
|
|
! decide whether or not to place a spin.
|
|
|
|
n_spin = 0
|
|
|
|
DO j = 0, nw-1
|
|
DO k = 0, nw-1
|
|
DO l = 0, nd-1
|
|
IF ( mod(j+k+l,2) .EQ. 0 ) THEN ! This takes care of the fcc structure.
|
|
IF ( ran(iseed) .LT. c ) THEN ! Takes care of concentration
|
|
n_spin = n_spin + 1
|
|
IF ( n_spin .GT. max_spins ) STOP ' Stopped because number of spin too large '
|
|
s(n_spin).x = j
|
|
s(n_spin).y = k
|
|
s(n_spin).z = l
|
|
|
|
! Give the spin an arbitrary direction
|
|
|
|
s(n_spin).theta = Pi * ran(iseed)
|
|
s(n_spin).phi = (Pi+Pi) * ran(iseed)
|
|
s(n_spin).dir(1) = sin(s(n_spin).theta) * cos(s(n_spin).phi)
|
|
s(n_spin).dir(2) = sin(s(n_spin).theta) * sin(s(n_spin).phi)
|
|
s(n_spin).dir(3) = cos(s(n_spin).theta)
|
|
|
|
END IF
|
|
END IF
|
|
|
|
|
|
END DO
|
|
END DO
|
|
END DO
|
|
|
|
! determine positions of the muons
|
|
|
|
n_site = 0
|
|
|
|
DO j = 1, nw-1, 2
|
|
DO k = 1, nw-1, 2
|
|
DO l = nd1, nd2, 2
|
|
|
|
IF ( ran(iseed) .LT. fraction ) THEN
|
|
n_site = n_site + 1
|
|
m(n_site).x = j
|
|
m(n_site).y = k
|
|
m(n_site).x = l
|
|
ns = 0
|
|
|
|
DO i = 1, n_spin
|
|
kw = j - s(i).x
|
|
IF ( kw .LT. -hw ) kw = kw + nw ! periodic boundary condition
|
|
IF ( kw .GT. hw ) kw = kw - nw ! periodic boundary condition
|
|
r(1) = dble(float(kw))
|
|
kw = k - s(i).y
|
|
IF ( kw .LT. -hw ) kw = kw + nw ! periodic boundary condition
|
|
IF ( kw .GT. hw ) kw = kw - nw ! periodic boundary condition
|
|
r(2) = dble(float(kw))
|
|
r(3) = dble(float(l-s(i).z)) ! NO periodic boundary condition
|
|
r_2 = sum( r * r )
|
|
|
|
IF ( r_2 .LE. radiussq ) THEN ! skip calculation if distance
|
|
! is too large
|
|
help = sqrt( r_2 )
|
|
r_3 = r_2 * help
|
|
r_5 = r_2 * r_3
|
|
ns = ns + 1
|
|
IF (ns .GT. max_nn) STOP ' Stopped because NS becomes too large '
|
|
m(n_site).s(ns) = i
|
|
m(n_site).r(1,ns) = r(1)
|
|
m(n_site).r(2,ns) = r(2)
|
|
m(n_site).r(3,ns) = r(3)
|
|
m(n_site).r_2(ns) = r_2
|
|
m(n_site).r_5(ns) = r_5
|
|
END IF
|
|
|
|
END DO
|
|
|
|
m(n_site).ns = ns
|
|
|
|
END IF
|
|
|
|
END DO
|
|
END DO
|
|
END DO
|
|
|
|
|
|
|
|
RETURN
|
|
END
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
! FIELDS calculates all internal fields at the muonsites and
|
|
! adds the external field
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
Subroutine fields( a, moment, b_ext, s, n_site, m)
|
|
|
|
implicit none
|
|
|
|
integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
|
|
& max_muons = 10000, & ! maximum number of muons
|
|
& max_nn = 500 ! maximum number of nearest neighbours
|
|
|
|
|
|
structure /spin/
|
|
integer*4 x,y,z
|
|
real*8 theta,phi,dir(3)
|
|
end structure
|
|
|
|
structure /muon/
|
|
integer*4 x,y,z, ns, s(max_nn)
|
|
real*8 dir(3), r(3,max_nn), r_2(max_nn), r_5(max_nn), omega(3)
|
|
end structure
|
|
|
|
real*8 Pi, Gyro, p_r, a, factor
|
|
real*8 b(3), b_ext(3), moment
|
|
integer*4 j, k, l, n_site
|
|
record /spin/ s(*)
|
|
record /muon/ m(*)
|
|
|
|
Pi = acos(-1D+00)
|
|
Gyro = (Pi+Pi) * 135.54 ! gyro-magnetic ratio of muon [tesla^-1 s^-1]
|
|
factor = 1D-07 * moment * 9.2740019D-24 / ( a*a*a * 0.125D-27 )
|
|
|
|
DO j = 1, n_site
|
|
b = 0
|
|
DO k = 1, m(j).ns
|
|
l = m(j).s(k)
|
|
p_r = m(j).r(1,k) * s(l).dir(1) + &
|
|
& m(j).r(2,k) * s(l).dir(2) + &
|
|
& m(j).r(3,k) * s(l).dir(3)
|
|
b(1) = b(1) + (3.0D+00*p_r*m(j).r(1,k)-m(j).r_2(k)*s(l).dir(1))/m(j).r_5(k)
|
|
b(2) = b(2) + (3.0D+00*p_r*m(j).r(2,k)-m(j).r_2(k)*s(l).dir(2))/m(j).r_5(k)
|
|
b(3) = b(3) + (3.0D+00*p_r*m(j).r(3,k)-m(j).r_2(k)*s(l).dir(3))/m(j).r_5(k)
|
|
END DO
|
|
b = factor * b + b_ext
|
|
m(j).omega(1) = Gyro * b(1)
|
|
m(j).omega(2) = Gyro * b(2)
|
|
m(j).omega(3) = Gyro * b(3)
|
|
END DO
|
|
|
|
RETURN
|
|
END
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
! FLUCTUATION changes all directions of the spins with a random amount
|
|
! DTHETA and DPHI
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
Subroutine fluctuation( iseed, n_spin, s, dtheta, dphi, fraction )
|
|
|
|
implicit none
|
|
|
|
integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
|
|
& max_muons = 10000, & ! maximum number of muons
|
|
& max_nn = 500 ! maximum number of nearest neighbours
|
|
|
|
structure /spin/
|
|
integer*4 x,y,z
|
|
real*8 theta,phi,dir(3)
|
|
end structure
|
|
|
|
record /spin/ s(*)
|
|
real*8 dtheta, dphi, dt, dp, Pi, fraction
|
|
integer*4 iseed, n_spin, i_spin
|
|
|
|
IF ( fraction .LT. 1.0D-06 .OR. &
|
|
& ( dtheta .LT. 1.0D-06 .AND. dphi .LT. 1.0D-06 ) ) RETURN
|
|
|
|
Pi = acos( -1.0D+00 )
|
|
dt = dtheta * Pi / 180.0 ! amplitude of the fluctuation in theta
|
|
dp = dphi * Pi / 180.0 ! amplitude of the fluctuation in phi
|
|
|
|
DO i_spin = 1, n_spin
|
|
IF ( ran(iseed) .LT. fraction ) THEN
|
|
s(i_spin).theta = s(i_spin).theta + 2.0 * dt * (ran(iseed)-0.5)
|
|
s(i_spin).phi = s(i_spin).phi + 2.0 * dp * (ran(iseed)-0.5)
|
|
s(i_spin).dir(1) = sin(s(i_spin).theta) * cos(s(i_spin).phi)
|
|
s(i_spin).dir(2) = sin(s(i_spin).theta) * sin(s(i_spin).phi)
|
|
s(i_spin).dir(3) = cos(s(i_spin).theta)
|
|
END IF
|
|
END DO
|
|
|
|
RETURN
|
|
END
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
! MUONROTATION rotates all muons over the vector m.omega * step
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
Subroutine muonrotation( n_site, m, step )
|
|
|
|
implicit none
|
|
|
|
integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
|
|
& max_muons = 10000, & ! maximum number of muons
|
|
& max_nn = 500 ! maximum number of nearest neighbours
|
|
|
|
structure /muon/
|
|
integer*4 x,y,z, ns, s(max_nn)
|
|
real*8 dir(3), r(3,max_nn), r_2(max_nn), r_5(max_nn), omega(3)
|
|
end structure
|
|
|
|
record /muon/ m(*)
|
|
real*8 v(3), OM(3), step
|
|
integer*4 j, n_site
|
|
|
|
DO j = 1, n_site
|
|
OM(1) = m(j).omega(1)
|
|
OM(2) = m(j).omega(2)
|
|
OM(3) = m(j).omega(3)
|
|
OM = OM * step
|
|
v(1) = m(j).dir(1)
|
|
v(2) = m(j).dir(2)
|
|
v(3) = m(j).dir(3)
|
|
call rotation( v, OM )
|
|
m(j).dir(1) = v(1)
|
|
m(j).dir(2) = v(2)
|
|
m(j).dir(3) = v(3)
|
|
END DO
|
|
|
|
RETURN
|
|
END
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
! ROTATION rotates a vector V around the vector O
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|
|
|
|
Subroutine rotation( v, o )
|
|
|
|
implicit none
|
|
|
|
real*8 v(3), o(3), uo(3), r(3), o_abs, cc, ss
|
|
|
|
o_abs = sqrt( sum( o * o ) )
|
|
|
|
IF ( o_abs .GT. 1D-08 ) THEN
|
|
uo = o / o_abs
|
|
cc = cos( o_abs )
|
|
ss = sin( o_abs )
|
|
|
|
r(1) = ( cc+uo(1)*uo(1)*(1-cc) ) * v(1) + &
|
|
& ( -uo(3)*ss+uo(1)*uo(2)*(1-cc) ) * v(2) + &
|
|
& ( uo(2)*ss+uo(1)*uo(3)*(1-cc) ) * v(3)
|
|
|
|
r(2) = ( uo(3)*ss+uo(1)*uo(2)*(1-cc) ) * v(1) + &
|
|
& ( cc+uo(2)*uo(2)*(1-cc) ) * v(2) + &
|
|
& ( -uo(1)*ss+uo(2)*uo(3)*(1-cc) ) * v(3)
|
|
|
|
r(3) = ( -uo(2)*ss+uo(1)*uo(3)*(1-cc) ) * v(1) + &
|
|
& ( uo(1)*ss+uo(2)*uo(3)*(1-cc) ) * v(2) + &
|
|
& ( cc+uo(3)*uo(3)*(1-cc) ) * v(3)
|
|
|
|
v = r
|
|
END IF
|
|
|
|
RETURN
|
|
END
|