(load "ncad021.scm")

(define (methane)
  (clear-structure)
  (add-atom "C" #( 0.0  0.0  0.0))
  (add-atom "H" #( 1.0  1.0  0.0))
  (add-atom "H" #( 2.0 -1.0  0.0))
  (add-atom "H" #(-1.0  0.0  1.0))
  (add-atom "H" #(-1.0  0.0 -1.0))
  (add-bond 1 0 1)
  (add-bond 1 0 2)
  (add-bond 1 0 3)
  (add-bond 1 0 4))

(define (water)
  (add-atom "O" #( 0.0 0.0 0.0))
  (add-atom "H" #(-1.0 1.0 0.0))
  (add-atom "H" #( 1.0 1.0 0.0))
  (add-bond 0 0 1)
  (add-bond 0 0 2))

(define (show-structure)
  (let ((atom-summary
	 (lambda (a)
	   (list ((a 'element) 'name)
		 (a 'position)))))
    (printf "~s~%" (atom-summary (car atom-list)))))

(define (cook-it)
  (show-structure)
  ;; coarse energy minimization
  (set! emin-factor coarse-emin-factor)
  (emin-step)
  (show-structure)
  (emin-step)
  (show-structure)
  ;; fine energy minimization
  (set! emin-factor fine-emin-factor)
  (emin-step)
  (show-structure)
  (emin-step)
  (show-structure))

;; how does propane deform if two of the hydrogens repel each other?
;; here's how to find out

(define (repulsion a1 a2 f)
  (let* ((p1 (a1 'position))
	 (p2 (a2 'position))
	 (diff (vdiff p1 p2))
	 (d (vlen diff))
	 (fvec (vscale diff (/ f d))))
    (a1 'add-force fvec)
    (a2 'add-force (vscale fvec -1))))

(define (external-forces)
  (repulsion (list-ref atom-list 3)
	     (list-ref atom-list 9) 0.5))

(load-structure "propane")
(cook-it)
(save-structure-xyz "propane.xyz")

;; if you do this in MrEd, the GUI will pop up and show you what the
;; deformed propane molecule looks like. otherwise, you need to look
;; at propane.xyz using RasMol or some other viewer.