Micro and Nano Mechanics Group

Create dislocations for ParaDiS input file

Sylvie Aubry (10/29/07)

     implicit none
     integer nmax
     parameter (nmax = 20)
     integer ndis,n
     integer i,j,t,k,m,inc,m0,m1,method
     integer tmin(nmax),tmax(nmax),tstep(nmax)
     integer length,fp,maxstep
     integer ncount,ncount_int(nmax+1),ncount_tmp(nmax+1)
     integer seed(nmax)
     double precision bx(nmax),by(nmax),bz(nmax)
     double precision nx(nmax),ny(nmax),nz(nmax)
     double precision a(nmax),b(nmax),c(nmax)
     double precision x0(nmax),y0(nmax),z0(nmax)
     double precision x,y,z,tmp
     double precision rc,rTol
     double precision xmin,ymin,zmin,xmax,ymax,zmax
     double precision bnorm,nnorm
     double precision burgMag,L,Lmin,Lmax,minSeg,maxSeg
     character*50 name,name2,type

Read in some parameters

L : Length in micrometer 

     read(*,*) L

maxstep for the simulation 

     read(*,*) maxstep

L = (read-in L) micrometer = 10^-6 m. 

     burgMag = 2.725d-10
     L = L*1d-6/burgMag 
     Lmin = -L/2.d0 + 300.d0
     Lmax =  L/2.d0
     write(*,*)
     write(*,*) 'L=',L
     write(*,*) 'Lmin=',Lmin
     write(*,*) 'Lmax=',Lmax
     write(*,*) 'Total number of steps=',maxstep

number of points on each dislocation line 

     inc = 50
     write(*,5) 'Number of points on the curve=',inc
     write(*,*)

File name 

     read(*,*) name

ndis = how many dislocations 

     read(*,*) ndis
     if (ndis.le.0.and.ndis.gt.nmax) then
        write(*,*) 'Wrong number of dislocations.'
        stop
     else
        write(*,*) 'Reading ',ndis,' dislocations.'
     endif

Read Burgers vectors and glide planes 

     do i=1,ndis
        write(*,*)
        write(*,30) 'For dislocation ',i,'/',ndis

Burgers vector 

        read(*,*) bx(i),by(i),bz(i)

Glide planes 

        read(*,*) nx(i),ny(i),nz(i)

Normalize Burgers vector/Glide plane in units of burgMag = a sqrt(3)/2. 

        do j=1,3
           bnorm = sqrt(bx(i)**2 + by(i)**2 + bz(i)**2)
           bx(i)= bx(i)/bnorm
           by(i)= by(i)/bnorm
           bz(i)= bz(i)/bnorm

           nnorm = sqrt(nx(i)**2 + ny(i)**2 + nz(i)**2)
           nx(i)= nx(i)/nnorm
           ny(i)= ny(i)/nnorm
           nz(i)= nz(i)/nnorm
        enddo

Choose how to get dislocation line 

        read(*,*) method
        if (method.eq.0) then

Define dislocation line automatically according to type and seed

Type of the dislocation : screw, edge, mixed 

           read(*,*) type
           write(*,*) 'Type of dislocation: ',type(1:length(type))
           read(*,*) seed(i)
   
           call get_abc(a(i),b(i),c(i),x0(i),y0(i),z0(i),
    &        bx(i),by(i),bz(i),nx(i),ny(i),nz(i),
    &           type,Lmin,Lmax,seed(i))
           
        else
           read(*,*) a(i),b(i),c(i),x0(i),y0(i),z0(i)
        endif

        write(*,20) 'bx=',bx(i),' by=',by(i),' bz=',bz(i)
        write(*,20) 'nx=',nx(i),' ny=',ny(i),' c=',nz(i)

        write(*,20) 'a=',a(i),' b=',b(i),' c=',c(i)
        write(*,20) 'x0=',x0(i),' y0=',y0(i),' z0=',z0(i)

        call get_bounds(xmin,xmax,Lmin,Lmax,a(i),x0(i))
        call get_bounds(ymin,ymax,Lmin,Lmax,b(i),y0(i))
        call get_bounds(zmin,zmax,Lmin,Lmax,c(i),z0(i))



        tmin(i) = max(xmin,ymin,zmin)
        tmax(i) = min(xmax,ymax,zmax)

        if (tmin(i).gt.tmax(i)) then
           tmp = tmin(i)
           tmin(i) = tmax(i)
           tmax(i) = tmp
        endif

inc in the final number of points on the dislocation line. 

        tstep(i) = (tmax(i)-tmin(i))/inc
        if (tstep(i).eq.0) then
           write(*,*) 'Found a number of point on the ',
    &                  'curve equals to zero - ABORT'
           stop
        endif
        write(*,*) 'bounds=',tmin(i),tmax(i)
     enddo

Write outputs files

write .cn file 

     name2 = 'Outputs/'//name(1:length(name))//'_results'
     open(unit=20,file='../../Runs/'//name(1:length(name))//'.cn')
     write(20,400) '### ParaDis input file (create_dislocation.f) ###'
     write(20,400) '### S.A. (10/29/07)'
     write(20,*)
     write(20,400) '#Directory to write output files'
     write(20,*) 'dirname =',name2(1:length(name2))
     write(20,*) 
     write(20,400)  '#Input files for PBC image stresses'
     write(20,400)  'Rijmfile = "Inputs/Rijm.cube.out"'
     write(20,400)  'RijmPBCfile = "Inputs/RijmPBC.cube.out"'
     write(20,*)
     write(20,400)  '#Total simulation step'
     write(20,*)  'maxstep=',maxstep
     write(20,*)
     write(20,400)  '#The total number of CPUs'
     write(20,400)  'numXdoms = 1'
     write(20,400)  'numYdoms = 1'
     write(20,400)  'numZdoms = 1'
     write(20,*)
     write(20,400)  '#Cells for dislocation grouping (>=3)'
     write(20,400)  'numXcells = 3'
     write(20,400)  'numYcells = 3'
     write(20,400)  'numZcells = 3'
     write(20,*)

     write(20,400)  '#Fast multipole method specs.'
     write(20,*)  'fmEnabled = 0  #disable fast multipole'
     write(20,*)  'fmMPOrder = 2'
     write(20,*)  'fmTaylorOrder = 4'
     write(20,*)  'fmCorrectionTbl = "Inputs/fm-ctab.m2.t4.dat"'
     write(20,*)

     write(20,400)  'timestepIntegrator = "backward-euler"'
     write(20,*)
     write(20,400)  '#Simulation box size'
     write(20,*)  'xBoundMin =',Lmin
     write(20,*)  'xBoundMax =',Lmax
     write(20,*)  'yBoundMin =',Lmin
     write(20,*)  'yBoundMax =',Lmax
     write(20,*)  'zBoundMin =',Lmin
     write(20,*)  'zBoundMax =',Lmax
     write(20,*)

     write(20,400)  '#Boundary conditions'
     write(20,*)  'xBoundType = 0'
     write(20,*)  'yBoundType = 0'
     write(20,*)  'zBoundType = 0'
     write(20,*)

     write(20,400)  '#Mobility law function'
     write(20,*)  'mobilityLaw = "BCC_0"'
     write(20,*)  'MobScrew =', 10.
     write(20,*)  'MobEdge  =', 10.
     write(20,*)  'MobClimb =', 1.e-4

max/min Seg = max/min length of dislocation segment. current dislocation segment length is inc 

     write(20,*)
     minSeg = 4*inc
     maxSeg = 2*minSeg
     write(20,400)  '#Discretization'
     write(20,*)  'maxSeg =',maxSeg
     write(20,*)  'minSeg =',minSeg
     write(20,*)

     write(20,400)  '#Maximum nodal displacement at each time step'
     write(20,400)  'rmax = 100'
     write(20,*)
     write(20,400)  '#Error tolerance in determining time step'

It should be such that rntol < 0.5 rc = 0.5 a. 

     rc = minSeg/sqrt(3.d0)*0.5
     rTol = 0.45*0.5*rc
     write(20,*)  'rTol =',rTol    
     write(20,*)
     write(20,400)  '#Maximum time step'
     write(20,*)  'maxDT = 1e+07'
     write(20,*)
     write(20,*)  '#Core cut-off radius'
     write(20,*)  'rc =',rc
     write(20,*)

     write(20,400)  '#Turn on elastic interaction'
     write(20,*)  'elasticinteraction = 1'
     write(20,*)

     write(20,400)  '#Applied stress in Pa (xx,yy,zz,yz,zx,xy)'
     write(20,*)  'appliedStress = [ 0 0 0 0 0 0 ]'
     write(20,*)

     write(20,400)  '#Loading'
     write(20,*)  'loadType = 1' 
     write(20,*)  'eRate = 1'
     write(20,*)  'edotdir = [1 0 0]'
     write(20,*)

     write(20,400)  '#Save files'
     write(20,*)  'savecn = 1'
     write(20,*)  'savecnfreq = 100'
     write(20,*)

     write(20,400)  '#Povray output'
     write(20,*)  'povray = 1'
     write(20,*)  'povrayfreq = 100'
     write(20,*)  'povraydt = -1.000000e+00'
     write(20,*)  'povraytime = 0.000000e+00'
     write(20,*)  'povraycounter = 0'
     write(20,*)  'saveprop = 1'
     write(20,*)  'savepropfreq = 1'
     write(20,*)
     write(20,*)  'preserveOldTags = 1'
     close(20)

Write .data file 

     open(unit=20,file='../../Runs/'//name(1:length(name))//'.data')
     write(20,400) '### ParaDis nodal file (create_dislocation.f)'
     write(20,400) '### S.A. (10/29/07)###'
     write(20,400)

     write(20,400) '#File version number'
     write(20,400) '2'
     write(20,400)

     write(20,400) '#Number of data file segments'
     write(20,400) '1'
     write(20,400) '#Minimum coordinate values (x, y, z)'
     write(20,*) -L/2,-L/2,-L/2
     write(20,400) '#Maximum coordinate values (x, y, z)'
     write(20,*) L/2,L/2,L/2
     write(20,*)

     write(20,400) '#Node count'
     ncount = 0
     do i=1,ndis+1
        ncount_int(i) = 0
     enddo

     do i=1,ndis
        do t=tmin(i),tmax(i),tstep(i)
           ncount = ncount + 1
           ncount_int(i+1) = ncount_int(i+1) + 1 
        enddo
     enddo

     write(20,*) ncount
     write(20,400) '#Domain geometry (x, y, z)'
     write(20,400) '1 1 1'
     write(20,*) 
     write(20,400) '#Domain boundaries (x, y, z)'
     write(20,*) -L/2
     write(20,*) -L/2
     write(20,*) -L/2
     write(20,*)  L/2
     write(20,*)  L/2
     write(20,*)  L/2
     write(20,*)

     write(20,400) '#node_tag, x, y, z, num_arms, constraint'
     write(20,400) '#arm_tag, burgx, burgy, burgz, nx, ny, nz'
     write(20,400) '#length in unit of burgMag'
     m = 0
     m0 = 0
     m1 = 0
     do i=1,ndis

first 

        m0 = m
        m1 = m1 + ncount_int(i)
        t=tmin(i)
        x = x0(i) + a(i)*t
        y = y0(i) + b(i)*t
        z = z0(i) + c(i)*t        
        write(20,250) '0,',m,' ',x,y,z,2,0
        write(20,300) '     0,',m0+ncount_int(i+1)-1,
    &        bx(i),by(i),bz(i)
        write(20,350) '       ',nx(i),ny(i),nz(i)
        write(20,300) '     0,',m+1,-bx(i),-by(i),-bz(i)
        write(20,350) '       ',nx(i),ny(i),nz(i)
        m = m + 1

all 

        do t=tmin(i)+tstep(i),tmax(i)-tstep(i),tstep(i)
           x = x0(i) + a(i)*t
           y = y0(i) + b(i)*t
           z = z0(i) + c(i)*t        
           write(20,250) '0,',m,' ',x,y,z,2,0
           write(20,300) '    0,',m-1,bx(i),by(i),bz(i)
           write(20,350) '      ',nx(i),ny(i),nz(i)
           write(20,300) '    0,',m+1,-bx(i),-by(i),-bz(i)
           write(20,350) '      ',nx(i),ny(i),nz(i)
           m = m + 1
        enddo

last node 

        t = tmax(i)
        x = x0(i) + a(i)*t
        y = y0(i) + b(i)*t
        z = z0(i) + c(i)*t        
        write(20,250) '0,',m,' ',x,y,z,2,0
        write(20,300) '    0,',m-1,bx(i),by(i),bz(i)
        write(20,350) '      ',nx(i),ny(i),nz(i)
        write(20,300) '    0,',m1,-bx(i),-by(i),-bz(i)
        write(20,350) '      ',nx(i),ny(i),nz(i)
        m = m + 1
     enddo

     close(20)

5    format(a,i3)
10   format(a,f9.2,a)
20   format(a,f9.2,a,f9.2,a,f9.2)
30   format(a,i2,a,i3)
100  format(3f15.3,i2,a)
200  format(2i5,6f10.3,a)

250  format(a,i4,a,3f12.2,2i2)
300  format(a,i4,3f12.2)
350  format(a,3f12.2)
400  format(a)

     end

     double precision function min(a,b,c)
     implicit none
     double precision a,b,c,min0     
     if (a.lt.b) then
        min0 = a
     else
        min0 = b
     endif
     if (min0.lt.c) then
        min = min0
     else
        min = c
     endif
     return min
     end

     double precision function max(a,b,c)
     implicit none
     double precision a,b,c,max0    
     if (a.gt.b) then
        max0 = a
     else
        max0 = b
     endif
     if (max0.gt.c) then
        max = max0
     else
        max = c
     endif       
     end   



     integer function length(str)
     character*(*) str
     n = len(str)
10   if (n.gt.0) then
       if (str(n:n).eq.' ') then
         n = n - 1
         goto 10
       endif
     endif
     length = n
     return
     end



     subroutine random_gen(ranval,min,max,seed)
     integer seed
     real*4 ranval_tmp
     double precision min,max,ranval      
     ranval_tmp = rand(seed)
     ranval = ranval_tmp*(max-min) + min
     end



    subroutine get_abc(a,b,c,x0,y0,z0,bx,by,bz,nx,ny,nz,
    &     type,Lmin,Lmax,seed)
     implicit none
     integer length
     integer seed,seedX0,seedY0,seedZ0,seedalpha
     double precision a,b,c
     double precision x0,y0,z0,bx,by,bz,nx,ny,nz,Lmin,Lmax
     double precision alpha
     character*20 type

     seedX0 = seed + 19
     seedY0 = seed + 395
     seedZ0 = seed + 34
     seedalpha = seed

     call random_gen(x0,Lmin/5d0,Lmax/5d0,seedX0)
     call random_gen(y0,Lmin/5d0,Lmax/5d0,seedY0)
     call random_gen(z0,Lmin/5d0,Lmax/5d0,seedZ0)
     call random_gen(alpha,-100.d0,100.d0,seedalpha)

screw : dislocation line is parallel to the Burgers vector. 

     if (type(1:length(type)).eq.'screw') then
        a = alpha * bx
        b = alpha * by
        c = alpha * bz
     endif

edge : dislocation line is normal to the Burgers vector and the glide plane 

     if (type(1:length(type)).eq.'edge') then
        a =  alpha * (by*nz - bz*ny)
        b = -alpha * (bx*nz - bz*nx)
        c =  alpha * (bx*ny - by*nx)
     endif

mixed : edge and screw component 

     if (type(1:length(type)).eq.'mixed') then
        a = alpha * bx + alpha * (by*nz - bz*ny)
        b = alpha * by - alpha * (bx*nz - bz*nx)
        c = alpha * bz + alpha * (bx*ny - by*nx)
     endif

     end

     subroutine get_bounds(xmin,xmax,Lmin,Lmax,a,x)
     implicit none
     integer i
     double precision xmin,xmax,Lmin,Lmax,a,x
        if (a.ne.0) then
           if (a.gt.0) then
              xmin = (Lmin-x)/a
              xmax = (Lmax-x)/a
           else
              xmin = (Lmax-x)/a
              xmax = (Lmin-x)/a 
           endif
        else
           xmin = -10d10
           xmax =  10d10
        endif
     end



     subroutine get_rand(seed,yy)
     implicit none
     integer i, number, old, seed, x, y
     double precision yy     
      number = 1000
      old = seed
      do i = 1, number
         x = MOD((57*old+1), 256)
         y = MOD((57*x+1), 256)
        old=y 
     enddo
     yy = dabs(y*1.d0/255.d0)
     end
Retrieved from "http://micro.stanford.edu/mediawiki/index.php?title=Create_Straight_Dislocations_for_ParaDiS_Input&oldid=649"
Member Login