! Program to calculate dipolar fields in spinglasses, ! their distribution and the depolarization of the muon ! ! Ge Nieuwenhuys, March, September, October 2005 ! ! October 12: periodic boundary conditions in y- z plane ! October 14: random number start randomly (based on clock) for ! batch calculations. ! October 14: output-file-names are automatically indexed. ! October 17: oversized the recordlength of the direct-accessfile for ! unknown, but apparently essential reasons. ! ! Spins are located on a fcc lattice ! ! nspin number of spins ! nsp number of spins asked ! d thickness ! a lattice constant ! ah half of lattice constant ! Use DFPORT ! library only needed for obtaining CPU-time Use DFLIB ! ! 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 dir(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 ! parameter( max_spins = 3000000, & ! maximum number of magnetic moments & gyro = 135.5, & ! gyromagnetic ratio of muon & twpi = 6.2831, & ! two times Pi & radius = 2.0, & ! maxinum distance [nm] for ! the dipole-field will be calculated & range = 10.0, & ! maximum absolute value of the field expected & mrange = 4000, & ! range of the integer histograms & nrange = 80 ) ! range of the normalized histograms ! character*10 dddd, tttt, zone character*4 file_index integer*4 dt(8), ifile, l_calc character*80 comment, calculation, line logical in_open, out_open, g_t_open, his_open, sgl, sgl_open integer*4 j,k,l,m,n, nsp, nspin, nat, id, ihist(3,-mrange:mrange) integer*4 iseed, maxfield, minfield, ihis, ibin, nd1, nd2, kd, ld, mh record /spin/ s(max_spins) real*8 d, concentration, c, dd(max_spins), w, depth1, depth2 real*8 px(max_spins),py(max_spins), pz(max_spins) real*8 b(3), factor, moment, help, r_3, r_5, r(3), p_r, sq_3, h(3) real*8 fraction, norm, aver_b(3), sigma_b(3), delta(3), anisotropy, b_ext(3) real*8 g_t(3,0:999), omega, b_abs, b_sq, ca_sq, his, radiussq real*4 runtime(2), start_time, end_time real*8 eb(3), emu(3), cc, ss, theta, phi ! Write(6,*) ' ' Write(6,*) ' ---------------------------------------------------------------------' Write(6,*) ' | Program field-calculation of muons due to random static spins |' Write(6,*) ' | Version of October 31, 2005 |' Write(6,*) ' | |' Write(6,*) ' | Input can also be read from an input file that should be named |' Write(6,*) ' | .inp and contain: |' Write(6,*) ' | |' Write(6,*) ' | ext. field(3) ,thickness, width, c, number_of_muons, |' Write(6,*) ' | lattice-constant [nm], magnetic moment [mu_B], |' Write(6,*) ' | initial-muon-direction(theta, phi)[degrees], |' Write(6,*) ' | (muon-positions from) depth1, (to) depth2 [nm], |' Write(6,*) ' | anisotropy [isotropic=1, planar <1, axial >1 |' Write(6,*) ' | (neg: ferromagnetic along the |' Write(6,*) ' | x - axis (anisotropy = -1.0) |' Write(6,*) ' | y - axis (anisotropy = -2.0) |' Write(6,*) ' | z - axis (anisotropy = -3.0) |' Write(6,*) ' | |' Write(6,*) ' | O R |' Write(6,*) ' | |' Write(6,*) ' | name of the .sgl file produced by |' Write(6,*) ' | MAKE SPINGLASS (starting on the first position), |' Write(6,*) ' | number_of_muons, |' Write(6,*) ' | initial-muon-direction(theta, phi)[degrees], |' Write(6,*) ' | (muon-positions from) depth1, (to) depth2 [nm], |' Write(6,*) ' | |' write(6,*) ' | Lines starting with ! (first position) are treated as comments. |' Write(6,*) ' | can be issued as a commandline parameter |' Write(6,*) ' ---------------------------------------------------------------------' ! ! files : ! open(9,file='\simulations\counter.his',status='old') read(9,*) ifile ! initialize outputfile counter ! ! write(6,*) ' iargc = ', iargc() 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=200 ) open(2,file=calculation(1:l_calc)//'.out',status='unknown',action='write') ! END IF ! inquire(1, opened = in_open ) inquire(2, opened = out_open ) ! ! initialization of randomumber generator ! iseed = 1234567 ! ! Get eventually other values from the iput file ! 111 IF (in_open) THEN ! ! Read everything from the input file, one line per calculation ! ifile = ifile + 1 ! increase outputfile number rewind(9) write(9,*) ifile ! store for next program 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') open(4,file=calculation(1:l_calc)//file_index//'.his',status='unknown',action='write') ! inquire(3, opened = g_t_open ) inquire(4, opened = his_open ) ! 112 read(1,'(a80)',end=999) line IF ( ( line(1:1) .GE. 'a' .AND. line(1:1) .LE. 'z' ) .OR. & & ( line(1:1) .GE. 'A' .AND. line(1:1) .LE. 'Z' ) ) THEN l = index( line, ' ') - 1 write(6,*) line(1:l) open(7,file=line(1:l)//'.sgl',status='old', & & access='direct',form='binary',recl=40,action='read',err=998) read(line(l+1:80),*,err=998,end=999) n_site, theta, phi, depth1, depth2 ELSE IF ( line(1:1) .EQ. '!' ) THEN write(2,'(a)') line GOTO 112 ELSE read(line,*,err=998,end=999) a, moment, b_ext, d, w, concentration, & & n_site, theta, phi, depth1, depth2, anisotropy END IF END IF ! ! Initialize randomnumber generator "randomly" ! call date_and_time( dddd, tttt, zone, dt ) DO i = 1, dt(8) ! number milliseconds on the clock dummy = rand(iseed) END DO ! ELSE ! ! put standard values in the case of on-line calculation ! for the lattice (4 nm), moment (2 uB), external field (0,0,0) and ! initial_muon_spin in y-direction ! ! a = 0.4 ! Assume 0.4 nanometer moment = 2.0 ! Assume 2 Bohrmagneton per spin b_ext = 0.0 ! No external field emu = 0.0 emu(2) = 1.0 ! initial muon direction along y-axis ! ! ! Ask size of the system ! 3 write(6,4) 4 format( ' What thickness [nm] (0=stop) ? '\) read(5,*,err=3) d IF ( d .LT. 0.0 ) GOTO 3 IF ( d .EQ. 0.0 ) THEN Write(6,*) ' ' STOP ' program terminated by operator' END IF ! 5 write(6,6) 6 format( ' What width [nm] ? '\) read(5,*,err=5) w IF ( w .LE. 0.0 ) GOTO 5 depth1 = 0.0 depth2 = w ! 7 write(6,8) 8 format( ' Which concentration [at.%] ? '\) read(5,*,err=7) concentration IF ( concentration .LE. 0.0 ) GOTO 7 ! ! Ask for the anisotropy. ! The random value of the direction cosin in the x-direction is multiplied ! by anisotropy before normalization ! 9 write(6,10) 10 format( ' The random value of the direction cosin in the x-direction'/ & & ' is multiplied by anisotropy before normalization'/ & & ' Anisotropy [isotrope == 1] ? '\) read(5,*,err=9) anisotropy ! 20 write(6,21) 21 format( ' Give value of the external field (x=perp to film,'/ & & ' y=initial_muon > '\) read(5,*,err=20) b_ext ! END IF ! end reading from input file / keyboard ! !---------------------------------------------------------------------------------------- ! Start calculation !---------------------------------------------------------------------------------------- call date_and_time( dddd, tttt, zone, dt ) ! ! If a spinglass has been simulated by MAKE SPINGLASS, then ! the .sgl file will be read, ELSE a random ! glass will be generated here. ! inquire(7, opened = sgl_open ) ! IF ( sgl_open ) THEN ! spin glass has been made read(7,rec=1) n,m,nspin,a,moment,T_glass read(7,rec=2) concentration,b_ext,steps_per_spin DO ispin = 1, nspin read(7,rec=ispin+2) s(ispin) END DO close(7) ! ELSE ! spin glass has NOT been made ! c = concentration / 100.0 ! ! Calculate the 'rounded' number of spins for a lattice n*m*m for ! the given concentration. ! n is the number of atoms (half unitcells) perpendicular ! to the layer (== x-direction). ! m is the size of the layer in the y- ad z-direction ! n = floor(2.0 * d / a ) + 2 m = floor(2.0 * w / a ) + 2 nat = m * m * n / 2 nspin = floor( nat * c ) ! IF (nspin .GE. max_spins ) THEN Write(6,*) ' ' Write(6,*) ' Too many spins: ', nspin IF ( out_open ) Write(2,*) ' Too many spins: ', nspin GOTO 111 END IF ! ! 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. ! nspin = 0 ! DO j = 0, n-1 DO k = 0, m-1 DO l = 0, m-1 IF ( mod(j+k+l,2) .EQ. 0 ) THEN ! This takes care of the fcc structure. IF ( ran(iseed) .LT. c ) THEN nspin = nspin + 1 s(nspin).x = j s(nspin).y = k s(nspin).z = l IF (anisotropy .GE. 0.0 ) THEN ! ! Give the spin an arbitrary direction ! DO i = 1, 3 h(i) = 2.0D+00 * ran(iseed) - 1.0D+00 END DO ! ! The anisotropy is taken care off by ! multiplying the direction cosine in ! the x-direction with ANOSOTROPY ! before normalizing the direction cosines. ! h(1) = anisotropy * h(1) norm = sum( h * h ) h = h / sqrt( norm ) ELSE h = 0.0 h(-int(anisotropy)) = 1.0 END IF s(nspin).dir = h ! END IF END IF END DO END DO END DO ! ! The sample has been grown now. ! Write(6,*) ' ' Write(6,*) 'The sample has been grown, calculation can start' Write(6,*) ' ' ! END IF ! Of reading ,calculation>.sgl or ! growing magnetic structure ! ! Now start the serious work. ! ! Use half of the lattice parameter as unit of length ! ah = a / 2.0 ! ! help for periodic boundary conditions ! mh = m / 2 ! ! the maximum distance squared in units of ah: ! radiussq = radius * radius / ( ah * ah ) ! ! Calculate factor to translate to the correct dimensions. ! ! factor is ( mu_o / 4 Pi ) * moment * mu_B / ( ah^3 ) ! -- ALL in MKS units -- ! so that the "field" can be calculated as ! 1/r^5 ( 3 * (s.dir *** r) * r - r^2 s.dir ), ! where s.dir is the unit vector to the direction of the magnetic moment, ! and *** stands for the dot-product. ! factor = 1D-07 * moment * 9.2740019D-24 / ( ah*ah*ah * 1D-27 ) ! ! see where the muons should go ! nd1 = floor( depth1 / ah ) nd2 = floor( depth2 / ah ) IF ( mod( nd1 , 2 ) .EQ. 0 ) nd1 = nd1 + 1 ! nd1 should be odd IF ( nd2 .LT. nd1 + 1 ) nd2 = nd1 + 1 ! ! calculate unit vector along the initial muon-spin direction ! emu(1) = sin( twpi * theta / 360.0 ) * cos( twpi * phi / 360.0) emu(2) = sin( twpi * theta / 360.0 ) * sin( twpi * phi / 360.0) emu(3) = cos( twpi * theta / 360.0 ) ! ! Ask the number of sites to calculated, about 10,000 is reasonable ! IF ( .NOT. in_open ) THEN ! read keyboard if no input file ! write(6,*) ' total number of muon-sites :', (m-1)*(m-1)*(nd2-nd1+1) / 8 write(6,*) ' ' 11 write(6,12) 12 format(' Give number of sites to be calculated > ' $) read(5,*,err=11) n_site ! END IF ! of reading keyboard ! fraction = dble( float(n_site) / float( (m-1)*(m-1)*(nd2-nd1+1)/8 )) ! ! make some space ! Write(6,*) ' ' Write(6,*) ' ' ! start_time = dtime(runtime) ! record the starttime ! ! Initialize the averages ! ib = 0 ! index of field calculation aver_b = 0 ! average of the field sigma_b = 0 ! average of the field squared hist = 0 ! histograms g_t = 0.0 ! initialize the line ! ! Assume the muon to be in the center of the fcc-cube ! DO j = nd1, nd2, 2 DO k = 1, m-1, 2 DO l = 1, m-1, 2 ! ! These do-loops run over all sites, which is probably too much (time consuming) ! Therefore select randomly sufficient (see above) fraction of ! the possible muon sites and calculate the dipolar field. ! IF ( ran(iseed) .LT. fraction ) THEN ! ! Calculate the field by running over all spins. ! In calculating the mutual distance, periodic boundaryconditions are applied ! in the y- and z-direction, but NOT in the x-direction, since that is supposed ! perpendicular to the film ! ! The field is only calculated when the distance is smaller then radius ! b = 0 ! DO i = 1, nspin r(1) = dble(float(j-s(i).x)) kd = k - s(i).y IF ( kd .LT. -mh ) kd = kd + m ! periodic boundary condition IF ( kd .GT. mh ) kd = kd - m ! periodic boundary condition r(2) = dble(float(kd)) ld = l - s(i).z IF ( ld .LT. -mh ) ld = ld + m ! periodic boundary condition IF ( ld .GT. mh ) ld = ld - m ! periodic boundary condition r(3) = dble(float(ld)) 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 h = s(i).dir p_r = sum( h * r ) b = b + ( 3.0D+00 * p_r * r - r_2 * h ) / r_5 END IF ! END DO ! ib = ib + 1 ! count the sites calculated. b = factor * b ! get correct dimensions aver_b = aver_b + b ! add the field to the averages sigma_b = sigma_b + b*b ! ! ! Count for histograms ! DO ih = 1, 3 ival = int( float(mrange) * b(ih) / range + 0.5D+00 ) IF ( abs(ival) .LE. mrange ) ihist(ih,ival) = ihist(ih,ival) + 1 END DO ! b = b + b_ext ! add external field b_sq = sum( b * b ) ! square of the field b_abs = sqrt( b_sq ) ! absolute value eb = b / b_abs ! unit vector omega = gyro * twpi * b_abs ! precession frequency ! ! Calculate the rotation of the muonspin for 1000 time-steps. ! The contribution to the asymmetry equals the components of the temporal ! muonspin, assuming the counters to be forward-backward, left-right ,and up-down, ! respectively. ! DO it = 0, 999 t = 1.0D-02 * dble(float(it)) cc = cos( omega * t ) ss = sin( omega * t ) ! g_t(1,it) = g_t(1,it) + & & ( cc+eb(1)*eb(1)*(1-cc)) * emu(1) + & & ( -eb(3)*ss+eb(1)*eb(2)*(1-cc)) * emu(2) + & & ( eb(2)*ss+eb(1)*eb(3)*(1-cc)) * emu(3) ! g_t(2,it) = g_t(2,it) + & & ( eb(3)*ss+eb(1)*eb(2)*(1-cc)) * emu(1) + & & ( cc+eb(2)*eb(2)*(1-cc)) * emu(2) + & & ( -eb(1)*ss+eb(2)*eb(3)*(1-cc)) * emu(3) ! g_t(3,it) = g_t(3,it) + & & ( -eb(2)*ss+eb(1)*eb(3)*(1-cc)) * emu(1) + & & ( eb(1)*ss+eb(2)*eb(3)*(1-cc)) * emu(2) + & & ( cc+eb(3)*eb(3)*(1-cc)) * emu(3) ! END DO ! IF ( mod(ib,1000) .EQ. 0 ) idummy = putc('#') ! END IF ! decision on fraction of muon sites END DO END DO END DO ! l, k, j loops ! ! Average over all calculaled sites. ! norm = dble( float(ib)) aver_b = aver_b / norm sigma_b = sqrt( (sigma_b - aver_b * aver_b ) / norm ) delta = gyro * sigma_b g_t = g_t / norm ! ! Renormalize histograms ! IF ( his_open ) THEN ! Should the histogram be calculated ?? Write(4,*) '-------------------------------------------------------' ! ! Check whether the maximum calculated field exceeds the range ! IF ( ihist(1,-mrange) .EQ. 0 .AND. ihist(1,mrange) .EQ. 0 .AND. & & ihist(2,-mrange) .EQ. 0 .AND. ihist(2,mrange) .EQ. 0 .AND. & & ihist(3,-mrange) .EQ. 0 .AND. ihist(3,mrange) .EQ. 0 ) THEN ! ! determine the range of fields found ! DO j = 1, 3 DO k = -mrange, mrange IF ( ihist(j, k) .GT. 0 ) maxfield = k IF ( ihist(j,-k) .GT. 0 ) minfield = -k END DO ! ! adjust binning of histogram and write values ! ibin = (maxfield - minfield) / nrange + 1 x = float(minfield) * range / float(mrange) step = range * float(ibin) / float(mrange) ! write(6,*) ' The field histogram vaues are: ' write(6,*) minfield, maxfield, ibin, x, step ! DO i = minfield, maxfield, ibin ihis = 0 DO k = 0, ibin-1 ihis = ihis + ihist(j,i+k) END DO his = float(ihis) / norm Write(4,'(2E16.6)') x, his x = x + step END DO ! Write(4,*) ' ' END DO ! ELSE Write(4,*) ' Fields exceed the maximum field for histogram calculation ' END IF END IF ! Histogram calculation ! end_time = dtime(runtime) ! write(6,*) ' ' write(2,100) comment(1:73),(dt(j),j=1,3),(dt(j),j=5,8) write(6,101) n*ah, m*ah write(6,301) nd1*ah, nd2*ah write(6,102) concentration write(6,103) anisotropy, int(-anisotropy) write(6,104) n_site write(6,304) theta, phi write(6,105) nspin write(6,106) aver_b write(6,107) sigma_b write(6,108) delta write(6,308) b_ext write(6,109) end_time - start_time ! ! Look whether data have to be written to file ! IF ( out_open ) THEN write(2,100) comment(1:73),(dt(j),j=1,3),(dt(j),j=5,8) write(2,101) n*ah, m*ah write(2,301) nd1*ah, nd2*ah write(2,102) concentration write(2,103) anisotropy, int(-anisotropy) write(2,104) n_site write(2,304) theta, phi write(2,105) nspin write(2,106) aver_b write(2,107) sigma_b write(2,108) delta write(2,308) b_ext write(2,109) end_time - start_time END IF ! 100 format(' '/' ',73('-')/' ',a73/' ',73('-')/ & & ' Calculation started ',i5,'-',i2,'-',i2, & & ' at ',2(i2,':'),i2,'.',i3/' ',73('-')/' ') 101 format(' sample = ', F6.1, ' nanometer thick, and ', F6.1, ' nanometer wide.') 102 format(' concentration = ', F12.1, ' at. %') 103 format(' anisotropy = ', E12.3,' (int) ',I2) 104 format(' number of muons = ', I12) 105 format(' number of spins = ', I12) 106 format(' average field = ', 3E12.3,' tesla') 107 format(' second moment = ', 3E12.3,' tesla') 108 format(' corres. delta = ', 3E12.3,' 1/microseconde') 109 format(' cpu_time = ', E12.3, ' seconds') 308 format(' ext. field = ', 3E12.3,' tesla') 301 format(' penetration from = ', F6.1,' to ',F6.1' nanometer.') 304 format(' initial muon spin, theta = ',f6.2,' phi = ', f6.2) ! ! Write G(t) if the file is open ! 500 IF ( g_t_open ) THEN ! DO k = 0, 999 write(3,'(3E20.6)') (g_t(id,k),id=1,3) ! output END DO ! END IF ! ! Go back to read new parameters ! GOTO 111 ! ! On error in input_file ! 998 Write(6,*) ' ' Write(6,*) ' There is an error in the input file. ' IF ( out_open ) Write(2,*) ' There is an error in the input file. ' ! 999 IF ( in_open ) close(1) IF ( out_open ) close(2) IF ( g_t_open ) close(3) IF ( his_open ) close(4) END ! ! End of program !------------------------------------------------------------------------------------------- ! ! Functions and Subroutines ! !------------------------------------------------------------------------------------------- real*8 FUNCTION length( v ) real*8 v(3) length = sqrt( sum( v * v ) ) RETURN END ! real*8 FUNCTION scalar_product( v, w ) real*8 v(3), w(3) scalar_product = sum( v * w ) RETURN END ! real*8 FUNCTION length_vector_product( v, w ) real*8 v(3), w(3), vp(3), length call vector_product( vp, v, w ) length_vector_product = length( vp ) RETURN END ! SUBROUTINE vector_product( vp, v, w ) real*8 v(3), w(3), vp(3) vp(1) = v(2) * w(3) - v(3) * w(2) vp(2) = v(3) * w(1) - v(1) * w(3) vp(3) = v(1) * w(2) - v(2) * w(1) RETURN END