mpb.scm.in

麻省理工的计算光子晶体的程序

; Copyright (C) 1999, 2000, 2001, 2002, Massachusetts Institute of Technology.
;
; 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, write to the Free Software
; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

; ****************************************************************

; Get the number of arguments to a function p.  However, some
; older versions of Guile (e.g. 1.2) do not support the 'arity
; property, and procedure-property just returns false.  In
; this case, we assume that the procedure returns 1 argument,
; as this is the most useful default for our purposes.  Sigh.

(define (procedure-num-args p) 
  (let ((arity (procedure-property p 'arity)))
    (if arity (car arity) 1)))

; ****************************************************************

(define-class material-type no-parent)

(define-class dielectric material-type
  (define-property epsilon no-default 'number))

(define (index n) (epsilon (* n n))) ; convenient substitute for epsilon

(define-class dielectric-anisotropic material-type
  (define-property epsilon-diag no-default 'vector3)
  (define-property epsilon-offdiag (vector3 0 0 0) 'cvector3)
  (define-property
    epsilon-offdiag-imag (vector3 0 0 0) 'vector3))

(define-class material-function material-type
  (define-property material-func no-default 'function
    (lambda (p) (= 1 (procedure-num-args p)))))

; use the solid geometry classes, variables, etcetera in libgeom:
; (one specifications file can include another specifications file)
(include "@LIBCTL_DIR@/utils/geom.scm")

; ****************************************************************

; eigensolver flags (grabbed from eigensolver.h by configure)

; first, we must define a function (pow2 n) to return 2^n:
(define (pow2 n) (if (<= n 0) 1 (* 2 (pow2 (- n 1)))))

@EIGS_FLAGS_SCM@ ; substituted by configure script

; ****************************************************************

; More input/output variables (besides those defined by libgeom, above).

(define-input-var k-points '() (make-list-type 'vector3))

(define-input-var num-bands 1 'integer)
(define-input-var tolerance 1.0e-7 'number positive?)
(define-input-var target-freq 0.0 'number (lambda (x) (>= x 0)))

(define-input-var mesh-size 3 'integer positive?)

(define-input-var epsilon-input-file "" 'string)

(define-input-var deterministic? false 'boolean)

; Eigensolver minutiae:
(define-input-var simple-preconditioner? false 'boolean)
(define-input-var eigensolver-flags EIGS_DEFAULT_FLAGS 'integer)
(define-input-var eigensolver-block-size -11 'integer)
(define-input-var eigensolver-nwork 3 'integer positive?)
(define-input-var eigensolver-davidson? false 'boolean)
(define-input-output-var eigensolver-flops 0 'number)

(define-output-var freqs (make-list-type 'number))
(define-output-var iterations 'integer)

(define-output-var parity 'string)

(define-input-var negative-epsilon-ok? false 'boolean)
(define (allow-negative-epsilon)
  (set! negative-epsilon-ok? true)
  (set! target-freq (/ 1 infinity)))

; ****************************************************************

; Definitions of external (C) functions:

; (init-params p true) initializes the geometry, etcetera, and does
; everything else that's needed to get ready for an eigenvalue
; calculation with parity p (see below).  This should be called
; after the input variables are changed.  If false is passed instead
; of true, fields from a previous run are retained, if possible, as a
; starting point for the eigensolver.
(define-external-function init-params true false
  no-return-value 'integer 'boolean)

; (set-parity p) changes the parity that is solved for by
; solve-kpoint, below.  p should be one of the following constants
; init-params should already have been called.  Be sure to call
; (randomize-fields) if you change the parity without calling
; init-params.
(define NO-PARITY 0)
(define EVEN-Z 1)
(define ODD-Z 2)
(define EVEN-Y 4)
(define ODD-Y 8)
(define TE EVEN-Z)
(define TM ODD-Z)
(define PREV-PARITY -1)
(define-external-function set-parity false false 
  no-return-value 'integer)
(define set-polarization set-parity) ; backwards compatibility

; (randomize-fields) initializes the fields to random values; should
; only be called after init-params.
(define-external-function randomize-fields false false no-return-value)

; (solve-kpoint kpoint) solves for the specified bands at the given k point.
; Requires that (init-params) has been called, and does not re-read the
; input variables, but does write the output vars.
(define-external-function solve-kpoint false true no-return-value 'vector3)

(define-external-function get-dfield false false no-return-value 'integer)
(define-external-function get-hfield false false no-return-value 'integer)
(define-external-function get-efield-from-dfield false false no-return-value)
(define-external-function get-epsilon false false no-return-value)
(define-external-function fix-field-phase false false no-return-value)
(define-external-function compute-field-energy false false 
	(make-list-type 'number))

(define-external-function get-epsilon-point false false 'number 'vector3)
(define-external-function get-epsilon-inverse-tensor-point false false 
  'cmatrix3x3 'vector3)
(define-external-function get-energy-point false false 'number 'vector3)
(define get-scalar-field-point get-energy-point)
(define-external-function get-bloch-field-point false false 'cvector3 'vector3)
(define-external-function get-field-point false false 'cvector3 'vector3)

(define-external-function compute-energy-in-dielectric false false
  'number 'number 'number)
(define-external-function compute-field-integral false false
  'cnumber 'function)
(define-external-function compute-energy-integral false false
  'number 'function)
(define-external-function compute-energy-in-object-list false false
  'number (make-list-type 'geometric-object))

(define-external-function output-field-to-file false false
  no-return-value 'integer 'string)  

(define-external-function mpi-is-master? false false 'boolean)
(define-external-function using-mpi? false false 'boolean)
(define-external-function mpi-num-procs false false 'integer)
(define-external-function mpi-proc-index false false 'integer)

(define-external-function get-kpoint-index false false 'integer)
(define-external-function set-kpoint-index false false
  no-return-value 'integer)

(define-external-function sqmatrix-size false false 'integer 'SCM)
(define-external-function sqmatrix-ref false false 'cnumber 
  'SCM 'integer 'integer)
(define-external-function get-eigenvectors false false 'SCM
  'integer 'integer)
(define-external-function set-eigenvectors false false no-return-value
  'SCM 'integer)
(define-external-function dot-eigenvectors false false 'SCM
  'SCM 'integer)
(define-external-function scale-eigenvector false false no-return-value
  'integer 'cnumber)
(define-external-function output-eigenvectors false false no-return-value
  'SCM 'string)
(define-external-function input-eigenvectors false false 'SCM 'string 'integer)
(define-external-function save-eigenvectors false false no-return-value
  'string)
(define-external-function load-eigenvectors false false no-return-value
  'string)

(define cur-field 'cur-field)
(define-external-function cur-field? false false 'boolean 'SCM)
(define-external-function rscalar-field-make false false 'SCM 'SCM)
(define-external-function cvector-field-make false false 'SCM 'SCM)
(define-external-function cvector-field-nonbloch! false false
  no-return-value 'SCM)
(define-external-function field-make false false 'SCM 'SCM)
(define-external-function fields-conform? false false 'boolean 'SCM 'SCM)
(define-external-function field-set! false false no-return-value 'SCM 'SCM)
(define (field-copy f) (let ((f' (field-make f))) (field-set! f' f) f'))
(define-external-function field-load false false no-return-value 'SCM)
(define-external-function field-mapL! false false no-return-value 'SCM 
  'function (make-list-type 'SCM))
(define (field-map! dest f . src) (apply field-mapL! (list dest f src)))
(define-external-function integrate-fieldL false false 'cnumber
  'function (make-list-type 'SCM))
(define (integrate-fields f . src) (apply integrate-fieldL (list f src)))
(define-external-function rscalar-field-get-point false false 'number 
  'SCM 'vector3)
(define-external-function cvector-field-get-point false false 'cvector3 
  'SCM 'vector3)
(define-external-function cvector-field-get-point-bloch false false 'cvector3 
  'SCM 'vector3)

; ****************************************************************

; Set print-ok? to whether or not we are the MPI master process.
; However, don't try this if we are running within gen-ctl-io,
; as it won't work.
(if (not (defined? 'output-source)) ; (a function defined by gen-ctl-io)
    (set! print-ok? (mpi-is-master?)))

(if (and (not (defined? 'output-source)) (using-mpi?))
    (set! interactive? false)) ; MPI doesn't support interactive mode

; ****************************************************************

; Utility function to display a comma-delimited list of data for the
; current k point, prefixed by data-name and the current parity.

(define (display-kpoint-data data-name data)
  (print parity data-name ":, " (get-kpoint-index))
  (map (lambda (d) (print ", " d)) data)
  (print "\n"))

; ****************************************************************

; Computing parities:

(define-external-function compute-zparities false false
  (make-list-type 'number))
(define-external-function compute-yparities false false
  (make-list-type 'number))

(define (display-zparities)
  (display-kpoint-data "zparity" (compute-zparities)))
(define (display-yparities)
  (display-kpoint-data "yparity" (compute-yparities)))

; ****************************************************************

; Computing group velocities:

(define-external-function compute-group-velocity-component false false
  (make-list-type 'number) 'vector3)

; Return a list of the group velocity vector3's, in the cartesian
; basis (and units of c):
(define (compute-group-velocities)
  (let ((vx (compute-group-velocity-component
	     (cartesian->reciprocal (vector3 1 0 0))))
	(vy (compute-group-velocity-component
	     (cartesian->reciprocal (vector3 0 1 0))))
	(vz (compute-group-velocity-component
	     (cartesian->reciprocal (vector3 0 0 1)))))
    (map (lambda (x y z) (vector3 x y z)) vx vy vz)))

; Define a band function to be passed to run, so that you can easily
; display the group velocities for each k-point.
(define (display-group-velocities)
  (display-kpoint-data "velocity" (compute-group-velocities)))

; ****************************************************************

; Add some predefined variables, for convenience:

(define vacuum (make dielectric (epsilon 1.0)))
(define air vacuum)

(define nothing (make material-type)) ; punches a "hole" through objects
				      ; to the default/background material

(define infinity 1.0e20) ; big number for infinite dimensions of objects

(set! default-material air)

; ****************************************************************

; The remainder of this file consists of Scheme convenience functions.

; ****************************************************************

; Function to convert a k-point k into an equivalent point in the
; first Brillouin zone (not necessarily the irreducible Brillouin zone):
(define (first-brillouin-zone k)
  (define (n k) (vector3-norm (reciprocal->cartesian k)))
  (define (try+ k v)
    (if (< (n (vector3+ k v)) (n k)) (try+ (vector3+ k v) v) k))
  (define (try k v) (try+ (try+ k v) (vector3- (vector3 0) v)))
  (define (try-all k)
    (try (try (try k (vector3 1 0 0)) (vector3 0 1 0)) (vector3 0 0 1)))
  (define (try-all&repeat k)
    (let ((knew (try-all k)))
      (if (< (n knew) (n k)) (try-all&repeat knew) k)))
  (let ((k0 (vector3- k (vector-map inexact->exact k))))
    (if (< (n k0) (n k)) (try-all&repeat k0) (try-all&repeat k))))

; functions to manipulate the fields; these are mainly convenient
; wrappers for the external functions defined previously.

(define (get-efield which-band)
  (get-dfield which-band)
  (get-efield-from-dfield))

(define-param filename-prefix "")

(define (output-field)
  (output-field-to-file -1 filename-prefix))

(define (output-field-x)
  (output-field-to-file 0 filename-prefix))

(define (output-field-y)
  (output-field-to-file 1 filename-prefix))

(define (output-field-z)
  (output-field-to-file 2 filename-prefix))

(define (output-epsilon)
  (get-epsilon)
  (output-field-to-file -1 filename-prefix))

(define (compute-energy-in-objects . objects)
  (compute-energy-in-object-list objects))

; ****************************************************************
; Functions to compute and output gaps, given the lists of frequencies
; computed at each k point.

; The band-range-data is a list if ((min . k-point) . (max . k-point))
; pairs, with each pair describing the frequency range of a band and
; the k-points where it achieves its maximum/minimum.  Here, we update
; this data with a new list of band frequencies, and return the new
; data.  If band-range-data is null or too short, the needed entries
; will be created.
(define (update-band-range-data band-range-data freqs k-point)
  (define (ubrd band-range-data freqs br-start)
    (if (null? freqs)
	(append (reverse br-start) band-range-data)
	(let ((br (if (null? band-range-data)
		      (cons (cons infinity -1) (cons (- infinity) -1))
		      (car band-range-data)))
	      (br-rest (if (null? band-range-data) '() (cdr band-range-data))))
	  (let ((newmin (if (< (car freqs) (caar br))
			    (cons (car freqs) k-point) (car br)))
		(newmax (if (> (car freqs) (cadr br))
			    (cons (car freqs) k-point) (cdr br))))
	    (ubrd br-rest (cdr freqs) 
		  (cons (cons newmin newmax) br-start))))))
  (ubrd band-range-data freqs '()))

; Output the band range data in a nice format:
(define (output-band-range-data br-data)
  (define (obr br i)
    (if (not (null? br))
	(begin
	  (print "Band " i " range: " (caaar br) " at " (cdaar br)
		 " to "  (cadar br) " at " (cddar br) "\n")
	  (obr (cdr br) (+ i 1)))))
  (obr br-data 1))

; Output any gaps in the given band ranges, and return a list
; of the gaps as a list of (percent freq-min freq-max) lists.
(define (output-gaps band-range-data)
  (define (ogaps br-cur br-rest i gaps)
    (if (null? br-rest)
	(reverse gaps)
	(if (>= (cadr br-cur) (caaar br-rest))
	    (ogaps (car br-rest) (cdr br-rest) (+ i 1) gaps)
	    (let ((gap-size (/ (* 200 (- (caaar br-rest) (cadr br-cur)))
			       (+ (caaar br-rest) (cadr br-cur)))))
	      (print "Gap from band " i " (" (cadr br-cur) ") to band "
		     (+ i 1) " (" (caaar br-rest) "), " gap-size "%\n")
	      (ogaps (car br-rest) (cdr br-rest) (+ i 1)
		     (cons (list gap-size (cadr br-cur) (caaar br-rest)) gaps))
	      ))))
  (if (null? band-range-data)
      '()
      (ogaps (car band-range-data) (cdr band-range-data) 1 '())))