(if (not (defined? 'atom-list))
    (load "ncad021.scm"))

;; Use add-atom and add-bond to make up some simple structures

(define (propane)
  (clear-structure)
  (add-atom "C" '#(-1.6192231896357 0.10913828710019 0.059089635045815))
  (add-atom "C" '#(0.057508530150221 -0.59251320772937 0.057957421682237))
  (add-atom "C" '#(1.5614160525151 0.42533696428484 -0.11545370543238))
  (add-atom "H" '#(-2.4286228415902 -0.78350252298058 -0.4367304383663))
  (add-atom "H" '#(-1.6275048913605 1.1771713105404 -0.64271318143432))
  (add-atom "H" '#(-1.9450207470767 0.33650493105546 1.2919456677319))
  (add-atom "H" '#(0.18595554441562 -1.1546247770247 1.1821974294535))
  (add-atom "H" '#(0.088500477802441 -1.3735084054718 -0.94323432063565))
  (add-atom "H" '#(1.7674216822209 1.0546257613843 0.98788251635771))
  (add-atom "H" '#(2.5279590432339 -0.41575445499099 -0.32991347593932))
  (add-atom "H" '#(1.4316103393248 1.2171261138323 -1.1110275484632))
  (add-bond 1 2 8)
  (add-bond 1 9 2)
  (add-bond 1 2 10)
  (add-bond 1 7 1)
  (add-bond 1 1 6)
  (add-bond 1 1 2)
  (add-bond 1 0 1)
  (add-bond 1 0 3)
  (add-bond 1 4 0)
  (add-bond 1 5 0))

(define (methane)
  (clear-structure)
  (add-atom "C" '#(0.24079083903194 0.067969485745148 0.31796190177503))
  (add-atom "H" '#(-0.68923293935545 0.12473891294566 -0.53483688483037))
  (add-atom "H" '#(0.80032995726744 -0.90112884074279 -0.2678684281233))
  (add-atom "H" '#(-0.17474477280629 -0.22445027528141 1.5127937938113))
  (add-atom "H" '#(0.97282391566095 1.1397731585413 0.35212142860335))
  (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 1 0 1)
  (add-bond 1 0 2))

(define (save-define file-name)
  (entering "save-structure")
  (if (not (null? file-name))
      (let ()
        (delete-file file-name)
        (let ((outf (open-output-file file-name)))
	  (fprintf outf "(define (glop)~%")
	  (map (make-lambda (a)
		 (fprintf outf
			 "  (add-atom ~a~a~a #(~a ~a ~a))~%"
			 #\"
			 (element-name (atm-element a))
			 #\"
                         (vector-ref (atm-position a) 0)
                         (vector-ref (atm-position a) 1)
                         (vector-ref (atm-position a) 2)))
	       atom-list)
	  (map (make-lambda (b)
		 (fprintf outf
			 "  (add-bond ~a ~a ~a))~%"
			 (bond-order b)
			 (bond-first b)
			 (bond-second b)))
	       bond-list)
	  (fprintf outf ")~%")
          (close-output-port outf)))))

(define (show-structure)
  (let ((atom-summary
	 (make-lambda (a)
	   (list (element-name (atm-element a))
		 (atm-position a)))))
    (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))

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

;; Let's set up propane, and then set the external force as a square
;; wave of repulsion/attraction between two hydrogens, and see what
;; the result looks like

(define use-propane true)
(define atom-1 '())
(define atom-2 '())

(if use-propane
    (let ()
      (propane)
      (set! atom-1 (list-ref atom-list 3))
      (set! atom-2 (list-ref atom-list 9)))
    (let ()
      (methane)
      (set! atom-1 (list-ref atom-list 1))
      (set! atom-2 (list-ref atom-list 2))))

(define delta-t 0.02)

(define (time-step)
  (compute-forces)
  (dolist (a atom-list)
	  (atm-add-pos a (vscale (atm-velocity a) delta-t))
	  (atm-add-velocity a
	     (vscale (atm-force a)
		     (/ delta-t (element-mass (atm-element a)))))))

(define (fmod x y)
  (- x (* y (floor (/ x y)))))

(define _t 0.0)

(define period 2.0)
(define time-limit 0.4)
;; To do interesting animations on MrEd, increase time-limit to 3 or 4

(define (external-forces)
  (repulsion atom-1 atom-2
	     (if (> (fmod _t period) (* 0.5 period)) 3 -3)))

;; if we're running MrEd, draw a picture, else print coords
(if (defined? 'update-display)
    (define (blab) (update-display true))
    (define (blab)
      (printf "~s~%" (atm-position (list-ref atom-list 3)))))
;      (printf "~s~%" (atm-velocity (list-ref atom-list 3)))
;      (printf "~s~%" (atm-force (list-ref atom-list 3)))))

(define (g)
  (set! _t 0.0)
  (if (not use-mred) (printf "~%"))
  (do ()
      ((>= _t time-limit))
    (time-step)
;    (printf "~s~%" _t)
    (blab)
    (set! _t (+ _t delta-t))))

(define (h)
  (load "hackv.scm")
  (g))