;;   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) false)

;; ******************************************************************
;; DATA STRUCTURES

(define atom-list '())
(define bond-list '())
(define term-list '())  ;; terms are zero-arg procedures

;; ******************************************************************
;; The Periodic Table, with hybridized species

(define (hybridize-carbon bond-count)
  (entering "hybridize-carbon")
  (cond ((= (bond-count-singles bond-count) 4) 'sp3)
	((= (bond-count-triples bond-count) 1) 'sp)
	((= (bond-count-doubles bond-count) 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 bond-count) 1) 'sp2)
	(else 'sp3)))

(define periodic-table
  (list
   (create-element "C"  1.94 19.925  4 'sp3 (func hybridize-carbon))
   (create-element "H"  1.5  1.674   1 false   false)
   (create-element "O"  1.74 26.565  2 'sp3 (func hybridize-oxygen))
   (create-element "N"  1.82 23.251  3 false   false)
   (create-element "F"  1.65 31.545  1 false   false)
   (create-element "Cl" 2.03 38.064  1 false   false)
   (create-element "Br" 2.18 131.038 1 false   false)
   (create-element "I"  2.32 210.709 1 false   false)
   (create-element "S"  2.11 53.087  2 false   false)
   (create-element "Si" 2.25 46.454  4 false   false)
   (create-element "LP" 1.2  0.0     1 false   false)
   (create-element "P"  2.11 51.464  3 false   false)))

(define (lookup-element name)
  (entering "lookup-element")
  (do ((L periodic-table (cdr L)))
      ((or (null? L)
	   (equal? name (element-name (car L))))
       (if (null? L) false (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
   (create-species "C"  'sp3 0.357 1)        ;; sp3
   (create-species "C"  'sp2 0.357 2)        ;; sp2 alkene
   (create-species "C"  'sp2 0.357 3)        ;; sp2 carbonyl
   (create-species "C"  'sp  0.357 4)        ;; sp acetylene
   (create-species "H"  false   0.382 5)        ;; hydrocarbon
   (create-species "O"  'sp2 0.406 6)        ;; C-O-[HC]
   (create-species "O"  'sp3 0.536 7)        ;; O carbonyl
   (create-species "N"  false   0.447 8)        ;; sp3
   (create-species "F"  false   0.634 11)       ;; flouride
   (create-species "Cl" false   1.95  12)       ;; chloride
   (create-species "Br" false   2.599 13)       ;; bromide
   (create-species "I"  false   3.444 14)       ;; iodide
   (create-species "S"  false   1.641 15)       ;; sulfide
   (create-species "Si" false   1.137 19)       ;; silane
   (create-species "LP" false   0.13  20)       ;; lone pair
   (create-species "H"  false   0.292 21)       ;; alcohol
   (create-species "C"  false   0.357 22)       ;; cyclopropane
   (create-species "P"  false   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 (species-name (car L)))
		(equal? hybrid (species-hybridization (car L)))))
       (if (null? L) false (car L)))))

;; given an atom and a bond-count, hybridize the atom correctly, and
;; assign it a new species

(define (hybridize-atom atm bond-count)
  (entering "hybridize-atom")
  (let* ((elem (atm-element atm))
	 (name (element-name elem))
	 (howto (element-how-to-hybridize elem))
	 (hybrid (if howto (funcall howto bond-count) false)))
    (atm-set-species atm
		     (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 (atm-position a))))
	(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 (atm-position a))
     (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 (atm-position a))
     (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) '())
        ((funcall test (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 true)
	(delete-bond m n)
	(set! bond-list (cons (create-bond order m n) bond-list)))))

(define (delete-bond m n)
  (entering "delete-bond")
  (set! need-to-resetup-terms true)
  (set! bond-list
	(remove-from-list-if-test
	 bond-list
	 (make-lambda (bond) (or (and (= m (bond-first bond))
				      (= n (bond-second bond)))
				 (and (= n (bond-first bond))
				      (= m (bond-second bond))))))))

(define (add-atom e pos)
  (entering "add-atom")
  (set! need-to-resetup-terms true)
  (set! atom-list
	(append
	 atom-list
	 (list
	  (create-atm (lookup-element e) pos)))))

(define (delete-atom n)
  (entering "delete-atom")
  (set! need-to-resetup-terms true)
  (set! atom-list (remove-from-list atom-list n))
  (set! bond-list (remove-from-list-if-test
		   bond-list
		   (make-lambda (bond) (or (= n (bond-first bond))
					   (= n (bond-second bond))))))
  (set! bond-list
	(map
	 (make-lambda (bond)
		      (if (> (bond-first bond) n)
			  (set! bond (create-bond (bond-order bond)
						(- (bond-first bond) 1)
						(bond-second bond))))
		      (if (> (bond-second bond 'second) n)
			  (set! bond (create-bond (bond-order bond)
						(bond-first bond)
						(- (bond-second bond) 1))))
	   bond)
	 bond-list)))

;; ******************************************************************
;; Loading and saving structures

(define (clear-structure)
  (set! atom-list '())
  (set! bond-list '())
  (set! need-to-resetup-terms true))

(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
	       (make-lambda (z)
		 (create-atm (lookup-element (car z))
			    (vector (cadr z) (caddr z) (cadddr z))))
	       (car s)))
	(set! bond-list
	      (map
	       (make-lambda (z)
		 (create-bond (car z) (cadr z) (caddr z)))
	       (cadr s))))))

(do-scheme
 (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 (make-lambda (z)
		   (list (element-name (atm-element z))
			 (vector-ref (atm-position z) 0)
			 (vector-ref (atm-position z) 1)
			 (vector-ref (atm-position z) 2)))
		 atom-list)
	    (map (make-lambda (z)
		   (list (bond-order z)
			 (bond-first z)
			 (bond-second z)))
		 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~%"
	    (element-name (atm-element L))
	    (vector-ref (atm-position L) 0)
	    (vector-ref (atm-position L) 1)
	    (- (vector-ref (atm-position L) 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 use-mred
    (let ()
      (load "gui.scm")
      (set! this-session (make-object session%))))