|
|
;; ncad021.scm - Top level stuff, data structures, atom/bond bookkeeping
;; Copyright (C) 1996,1997 Will Ware
;;
;; 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., 675 Mass Ave, Cambridge, MA 02139, USA.
;;
;; I can be reached via email at .
;; These things are redefined in GUI
(define (error-msg txt) '())
(define (warning-msg txt) '())
(define this-session '())
(define (force-rotate-mode) '())
;; Hook for defining forces other than chemical
(define external-forces #f)
(define debugging #f)
(if debugging
(define-macro dbgprintf
(lambda (x . y)
`(printf ,x ,@y)))
(define-macro dbgprintf
(lambda (x . y) '())))
(if debugging
(define-macro entering
(lambda (name)
`(printf "Entering ~s~%" ,name)))
(define-macro entering
(lambda (name) '())))
;; ******************************************************************
;; Handy macros stolen from Common Lisp
(define-macro dolist
(lambda (args . body)
`(do ((local-list ,(cadr args) (cdr local-list))
(,(car args) '()))
((null? local-list))
(set! ,(car args) (car local-list))
,@body)))
(define-macro dotimes
(lambda (args . body)
`(do ((,(car args) 0 (+ 1 ,(car args)))
(iteration-limit ,(cadr args)))
((>= ,(car args) iteration-limit))
,@body)))
;; ******************************************************************
;; DATA STRUCTURES
(define atom-list '())
(define bond-list '())
(define term-list '()) ;; terms are zero-arg procedures
(define (make-atom elem pos)
(entering "make-atom")
(let ((force #(0.0 0.0 0.0))
(species (lookup-species (elem 'name)
(elem 'initial-hybridization))))
(lambda (x . y)
(case x
('element elem)
('set-element (set! elem (car y)))
('species species)
('set-species (set! species (car y)))
('position pos)
('add-pos (set! pos (vplus pos (car y))))
('set-pos (set! pos (car y)))
('force force)
('zero-force (set! force #(0.0 0.0 0.0)))
('add-force (set! force (vplus force (car y))))
(else
(printf "Atom trouble, args: ~s, ~s~%" x y))))))
(define (make-bond order first second)
(entering "make-bond")
(lambda (x)
(case x
('order order)
('first first)
('second second))))
(define (make-element name rvdw mass bonds init-hybrid ch-hybrid)
(entering "make-element")
(lambda (x)
(case x
('name name)
('rvdw rvdw)
('mass mass)
('total-bonds bonds)
('initial-hybridization init-hybrid)
('how-to-hybridize ch-hybrid))))
(define (make-species name hybrid evdw mm2)
(entering "make-species")
(lambda (x)
(case x
('name name)
('hybridization hybrid)
('evdw evdw)
('mm2-index mm2))))
(define (make-bond-count n)
(entering "make-bond-count")
(let ((singles 0)
(doubles 0)
(triples 0))
(dolist (bond bond-list)
(if (or (= n (bond 'first)) (= n (bond 'second)))
(case (bond 'order)
(1 (set! singles (+ singles 1)))
(2 (set! doubles (+ doubles 1)))
(else (set! triples (+ triples 1))))))
(lambda (x)
(case x
('singles singles)
('doubles doubles)
('triples triples)
('total-bonds (+ singles (* 2 doubles) (* 3 triples)))))))
;; ******************************************************************
;; The Periodic Table, with hybridized species
(define (hybridize-carbon bond-count)
(entering "hybridize-carbon")
(cond ((= (bond-count 'singles) 4) 'sp3)
((= (bond-count 'triples) 1) 'sp)
((= (bond-count 'doubles) 2) 'sp)
(else 'sp2)))
;; when ionized, oxygen can form a triple bond, but ignore that for now...
(define (hybridize-oxygen bond-count)
(entering "hybridize-oxygen")
(cond ((= (bond-count 'doubles) 1) 'sp2)
(else 'sp3)))
(define periodic-table
(list
(make-element "C" 1.94 19.925 4 'sp3 hybridize-carbon)
(make-element "H" 1.5 1.674 1 #f #f)
(make-element "O" 1.74 26.565 2 'sp3 hybridize-oxygen)
(make-element "N" 1.82 23.251 3 #f #f)
(make-element "F" 1.65 31.545 1 #f #f)
(make-element "Cl" 2.03 38.064 1 #f #f)
(make-element "Br" 2.18 131.038 1 #f #f)
(make-element "I" 2.32 210.709 1 #f #f)
(make-element "S" 2.11 53.087 2 #f #f)
(make-element "Si" 2.25 46.454 4 #f #f)
(make-element "LP" 1.2 0.0 1 #f #f)
(make-element "P" 2.11 51.464 3 #f #f)))
(define (lookup-element name)
(entering "lookup-element")
(do ((L periodic-table (cdr L)))
((or (null? L)
(equal? name ((car L) 'name)))
(if (null? L) #f (car L)))))
;; I've mostly succeeded in hiding MM2 as an implementation detail in
;; forces.scm, but there need to be a few concessions in this file.
(define mm2-table
(list
(make-species "C" 'sp3 0.357 1) ;; sp3
(make-species "C" 'sp2 0.357 2) ;; sp2 alkene
(make-species "C" 'sp2 0.357 3) ;; sp2 carbonyl
(make-species "C" 'sp 0.357 4) ;; sp acetylene
(make-species "H" #f 0.382 5) ;; hydrocarbon
(make-species "O" 'sp2 0.406 6) ;; C-O-[HC]
(make-species "O" 'sp3 0.536 7) ;; O carbonyl
(make-species "N" #f 0.447 8) ;; sp3
(make-species "F" #f 0.634 11) ;; flouride
(make-species "Cl" #f 1.95 12) ;; chloride
(make-species "Br" #f 2.599 13) ;; bromide
(make-species "I" #f 3.444 14) ;; iodide
(make-species "S" #f 1.641 15) ;; sulfide
(make-species "Si" #f 1.137 19) ;; silane
(make-species "LP" #f 0.13 20) ;; lone pair
(make-species "H" #f 0.292 21) ;; alcohol
(make-species "C" #f 0.357 22) ;; cyclopropane
(make-species "P" #f 1.641 25))) ;; phosphide
;; given an element name and a hybridization, look up the species
;; in mm2-table
(define (lookup-species name hybrid)
(entering "lookup-species")
(do ((L mm2-table (cdr L)))
((or (null? L)
(and (equal? name ((car L) 'name))
(equal? hybrid ((car L) 'hybridization))))
(if (null? L) #f (car L)))))
;; given an atom and a bond-count, hybridize the atom correctly, and
;; assign it a new species
(define (hybridize-atom atom bond-count)
(entering "hybridize-atom")
(let* ((elem (atom 'element))
(name (elem 'name))
(howto (elem 'how-to-hybridize))
(hybrid (if howto (howto bond-count) #f)))
(atom 'set-species
(lookup-species name hybrid))))
;; ******************************************************************
;; Rotations and Centering
(define (center-structure)
(entering "center-structure")
(if (> (length atom-list) 0)
(let ((cog (vector 0.0 0.0 0.0))) ;; center of gravity
(dolist (a atom-list)
(set! cog
(vplus cog (a 'position))))
(set! cog
(vscale cog (/ -1.0 (length atom-list))))
(dolist (a atom-list)
(a 'add-pos cog)))))
(define (rotate-all-atoms-x-axis theta)
(entering "rotate-all-atoms-x-axis")
(let ((ct (cos theta))
(st (sin theta))
(temp 0.0)
(tmpv #(0.0 0.0 0.0)))
(dolist
(a atom-list)
(set! tmpv (a 'position))
(set! temp (- (* ct (vector-ref tmpv 1))
(* st (vector-ref tmpv 2))))
(vector-set! tmpv 2
(+ (* ct (vector-ref tmpv 2))
(* st (vector-ref tmpv 1))))
(vector-set! tmpv 1 temp)
(a 'set-pos tmpv))))
(define (rotate-all-atoms-y-axis theta)
(entering "rotate-all-atoms-y-axis")
(let ((ct (cos theta))
(st (sin theta))
(temp 0.0)
(tmpv #(0.0 0.0 0.0)))
(dolist
(a atom-list)
(set! tmpv (a 'position))
(set! temp (+ (* ct (vector-ref tmpv 0))
(* st (vector-ref tmpv 2))))
(vector-set! tmpv 2
(- (* ct (vector-ref tmpv 2))
(* st (vector-ref tmpv 0))))
(vector-set! tmpv 0 temp)
(a 'set-pos tmpv))))
(define negligible-angle-sq 0.00001)
(define (rotate-structure xt yt)
(entering "rotate-structure")
(if (> (* xt xt) negligible-angle-sq)
(rotate-all-atoms-y-axis xt))
(if (> (* yt yt) negligible-angle-sq)
(rotate-all-atoms-x-axis yt)))
;; ******************************************************************
;; Add/delete/select atoms and bonds
(define (remove-from-list L n)
(cond ((null? L) '())
((= n 0) (cdr L))
(else (cons (car L)
(remove-from-list (cdr L) (- n 1))))))
(define (remove-from-list-if-test L test)
(cond ((null? L) '())
((apply test (list (car L)))
(remove-from-list-if-test (cdr L) test))
(else (cons (car L)
(remove-from-list-if-test (cdr L) test)))))
(define (add-bond order m n)
(entering "add-bond")
(if (not (= m n))
(let ()
(set! need-to-resetup-terms #t)
(delete-bond m n)
(set! bond-list (cons (make-bond order m n) bond-list)))))
(define (delete-bond m n)
(entering "delete-bond")
(set! need-to-resetup-terms #t)
(set! bond-list
(remove-from-list-if-test
bond-list
(lambda (bond) (or (and (= m (bond 'first))
(= n (bond 'second)))
(and (= n (bond 'first))
(= m (bond 'second))))))))
(define (add-atom e pos)
(entering "add-atom")
(set! need-to-resetup-terms #t)
(set! atom-list
(append
atom-list
(list
(make-atom (lookup-element e) pos)))))
(define (delete-atom n)
(entering "delete-atom")
(set! need-to-resetup-terms #t)
(set! atom-list (remove-from-list atom-list n))
(set! bond-list (remove-from-list-if-test
bond-list
(lambda (bond) (or (= n (bond 'first))
(= n (bond 'second))))))
(set! bond-list
(map
(lambda (bond)
(if (> (bond 'first) n)
(set! bond (make-bond (bond 'order)
(- (bond 'first) 1)
(bond 'second))))
(if (> (bond 'second) n)
(set! bond (make-bond (bond 'order)
(bond 'first)
(- (bond 'second) 1))))
bond)
bond-list)))
;; ******************************************************************
;; Loading and saving structures
(define (clear-structure)
(set! atom-list '())
(set! bond-list '())
(set! need-to-resetup-terms #t))
(define (load-structure file-name)
(entering "load-structure")
(clear-structure)
(if (not (null? file-name))
(let* ((inf (open-input-file file-name))
(s (read inf)))
(close-input-port inf)
(set! atom-list
(map
(lambda (z)
(make-atom (lookup-element (car z))
(vector (cadr z) (caddr z) (cadddr z))))
(car s)))
(set! bond-list
(map
(lambda (z)
(make-bond (car z) (cadr z) (caddr z)))
(cadr s))))))
(if (defined? 'pretty-print)
(define (print-to-file outf s)
(pretty-print s outf))
(define (print-to-file outf s)
(fprintf outf "~s~%" s)))
(define (save-structure file-name)
(entering "save-structure")
(if (not (null? file-name))
(let ()
(delete-file file-name)
(let ((outf (open-output-file file-name)))
(print-to-file
outf
(list
(map (lambda (z)
(list ((z 'element) 'name)
(vector-ref (z 'position) 0)
(vector-ref (z 'position) 1)
(vector-ref (z 'position) 2)))
atom-list)
(map (lambda (z)
(list (z 'order)
(z 'first)
(z 'second)))
bond-list)))
(close-output-port outf)))))
(define (save-structure-xyz file-name)
(entering "save-structure-xyz")
(if (not (null? file-name))
(let ()
(delete-file file-name)
(let ((outf (open-output-file file-name)))
(fprintf outf "~a~%Gray Goo~%"
(length atom-list))
(dolist
(L atom-list)
(fprintf
outf "~a ~a ~a ~a~%"
((L 'element) 'name)
(vector-ref (L 'position) 0)
(vector-ref (L 'position) 1)
(- (vector-ref (L 'position) 2))))
(close-output-port outf)))))
;; ******************************************************************
;; Fire up the GUI
(load "forces.scm")
;; If we're in MrEd, launch the GUI. If we're in MzScheme, don't bother
(if (defined? 'mred:startup)
(let ()
(load "gui.scm")
(set! this-session (make-object session%))))
|
