CCL Home Page
Up Directory CCL forces
;;   forces.scm - compute forces acting on atoms, based on MM2 energy terms
;;   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 .

;; ******************************************************************
;; Table of difference vectors and distances

;; At the beginning of compute-forces, we compute the difference vectors and
;; distances for all pairs of atoms. These would otherwise be repeatedly
;; calculated by various energy/force terms, so it's a big win to precalculate
;; all that stuff in a big batch and make it available for other functions.
;; For N atoms, we will need to compute N*(N-1)/2 vector/distance pairs.
;; But to keep my head from collapsing, I'll keep it simple and use N^2
;; things, the bookkeeping is too ugly otherwise.

(define diff-dist-table '())

(define (make-diff-dist m n)
  (let* ((ma (m 'position))
	 (na (n 'position))
	 (diff (vdiff ma na))
	 (dist (vlen diff)))
    (lambda (x)
      (case x
	('diff diff)
	('mdiff (vscale diff -1.0))
	('distance dist)
	(else (printf "(diff-dist ~s) ???~%" x))))))

(define (lookup-diff-dist m n)
  (vector-ref (vector-ref diff-dist-table n) m))

(define (build-diff-dist-table)
  (set! diff-dist-table (make-vector (length atom-list)))
  (dotimes
   (i (length atom-list))
   (let ((row (make-vector (length atom-list))))
     (dotimes
      (j (length atom-list))
      (vector-set! row j
		   (make-diff-dist (list-ref atom-list j)
				   (list-ref atom-list i))))
     (vector-set! diff-dist-table i row))))

;; ******************************************************************
;; Coefficients for MM2 energy/force terms

(define (lookup-helper compare-func result-func default-result the-list)
  (cond ((null? the-list) default-result)
	((compare-func (car the-list)) (result-func (car the-list)))
	(else (lookup-helper compare-func result-func default-result
			     (cdr the-list)))))

;; A length-coefficient entry is two MM2 atom types, followed by a
;; k_s value, followed by an r0 value.

(define length-coefficients
  '((1 5 460 1.113)
    (1 1 440 1.523)
    (2 2 960 1.337)
    (4 4 1560 1.212)
    (1 6 536 1.402)
    (1 8 510 1.438)
    (3 7 1080 1.208)
    (1 11 510 1.392)
    (1 12 323 1.795)
    (1 13 230 1.949)
    (1 14 220 2.149)
    (8 20 610 0.6)
    (8 8 560 1.381)
    (6 20 460 0.6)
    (6 21 460 0.942)
    (6 6 781 1.470)
    (1 19 297 1.880)
    (1 25 291 1.856)
    (1 15 321.3 1.815)
    (19 19 185 2.332)
    (22 22 440 1.501)))

(define (lookup-length-coeffs m n)
  (lookup-helper
   (lambda (x) (or (and (= m (car x))
			(= n (cadr x)))
		   (and (= n (car x))
			(= m (cadr x)))))
   (lambda (x) (cddr x))
   '(400 1.3)
   length-coefficients))

;; An angle-coefficient entry is three MM2 atoms types, followed by
;; a k_th value, followed by a th0 value.

(define angle-coefficients
  '((1 1 1 0.45 1.911)
    (1 1 5 0.36 1.909)
    (5 1 5 0.32 1.909)
    (1 1 11 0.65 1.911)
    (11 1 11 1.07 1.869)
    (1 2 1 0.45 2.046)
    (2 1 2 0.45 1.911)
    (1 2 2 0.55 2.119)
    (2 2 5 0.36 2.094)
    (2 2 2 0.43 2.094)
    (1 4 4 0.2 3.142)
    (1 3 7 0.46 2.138)
    (1 6 1 0.77 1.864)
    (1 8 1 0.63 1.88)
    (1 1 6 0.57 1.911)
    (1 6 20 0.35 1.835)
    (1 8 20 0.5 1.906)
    (20 6 20 0.24 2.286)
    (19 19 19 0.35 1.943)
    (19 1 19 0.4 2.016)
    (1 19 1 0.48 1.934)
    (12 1 12 1.08 1.95)
    (1 1 15 0.55 1.902)
    (1 15 1 0.72 1.902)
    ;; I think I made up these last few entries, they should be
    ;; considered suspect. Really, I should pick up the more complete
    ;; tables from Jay Ponder's ftp site.
    (4 4 5 0.4 3.142)
    (7 2 1 0.4 2.0944)
    (7 2 2 0.4 2.0944)
    (5 8 5 1.1 1.85)
    (1 2 6 0.4 2.0944)
    (1 2 5 0.4 2.0944)
    (5 2 6 0.4 2.0944)))

(define (lookup-angle-coeffs m n p)
  (lookup-helper
   (lambda (x) (or (and (= m (car x))
			(= n (cadr x))
			(= p (caddr x)))
		   (and (= p (car x))
			(= n (cadr x))
			(= m (caddr x)))))
   (lambda (x) (cdddr x))
   '(0.4 1.9)
   angle-coefficients))

;; An torsion-coefficient entry is four MM2 atoms types, followed by
;; three values for v1, v2, and v3 respectively.

(define torsion-coefficients
  '((1 1 1 1 1.39 1.88 0.65)
    (1 1 1 5 0.0 0.0 1.85)
    (5 1 1 5 0.0 0.0 1.65)
    (1 1 1 11 0.0 -0.6 6.46)
    (1 2 2 1 -0.69 69.47 0.0)
    (2 2 2 2 -6.46 55.58 0.0)
    (5 2 2 5 0.0 104.21 0.0)))

(define (lookup-torsion-coeffs m n p q)
  (lookup-helper
   (lambda (x) (or (and (= m (car x))
			(= n (cadr x))
			(= p (caddr x))
			(= q (cadddr x)))
		   (and (= q (car x))
			(= p (cadr x))
			(= n (caddr x))
			(= m (cadddr x)))))
   (lambda (x) (cddddr x))
   '(0.0 0.0 0.0)
   torsion-coefficients))

;; ******************************************************************
;; Geometric fun with vectors

(define (vplus v1 v2)
  (vector (+ (vector-ref v1 0) (vector-ref v2 0))
	  (+ (vector-ref v1 1) (vector-ref v2 1))
	  (+ (vector-ref v1 2) (vector-ref v2 2))))

(define (vdiff v1 v2)
  (vector (- (vector-ref v1 0) (vector-ref v2 0))
	  (- (vector-ref v1 1) (vector-ref v2 1))
	  (- (vector-ref v1 2) (vector-ref v2 2))))

(define (dot-product v1 v2)
  (+ (* (vector-ref v1 0) (vector-ref v2 0))
     (* (vector-ref v1 1) (vector-ref v2 1))
     (* (vector-ref v1 2) (vector-ref v2 2))))

(define (cross-product v1 v2)
  (vector (- (* (vector-ref v1 1) (vector-ref v2 2))
	     (* (vector-ref v1 2) (vector-ref v2 1)))
	  (- (* (vector-ref v1 2) (vector-ref v2 0))
	     (* (vector-ref v1 0) (vector-ref v2 2)))
	  (- (* (vector-ref v1 0) (vector-ref v2 1))
	     (* (vector-ref v1 1) (vector-ref v2 0)))))

(define (safe-sqrt x)
  (real-part (sqrt x)))

(define (safe-acos z)
  (if (> z 1.0) (set! z 1.0))
  (if (< z -1.0) (set! z -1.0))
  (acos z))

(define (vlen v)
  (safe-sqrt (dot-product v v)))

(define (vscale v x)
  (vector (* (vector-ref v 0) x)
	  (* (vector-ref v 1) x)
	  (* (vector-ref v 2) x)))

;; find the component of v1 perpendicular to v2
(define (perpendicular-component v1 v2)
  (vdiff v1 (vscale v2 (/ (dot-product v1 v2) (dot-product v2 v2)))))

(define (v-angle dd1 dd2)
  (safe-acos
   (/ (dot-product (dd1 'diff) (dd2 'diff))
      (* (dd1 'distance) (dd2 'distance)))))

;; i-torsion uses an obsolete version of v-angle, don't use torsion forces
;; for a while til I fix it

(define (i-torsion v1 v2 v3 v4)
  (let* ((d12 (lookup-diff-dist v1 v2))
	 (d23 (lookup-diff-dist v2 v3))
	 (d43 (lookup-diff-dist v4 v3))
	 (p (d12 'diff))
	 (q (d23 'diff))
	 (r (d43 'diff)))
    (v-angle (perpendicular-component p q)
	     (perpendicular-component r q))))

;; ******************************************************************
;; Forces computed by energy terms

;; Functions with names like 'build-???-term' run during setup-terms, and
;; they build zero-argument functions which are stored in the term list.
;; setup-terms runs fairly infrequently (only when there are changes in
;; the structure: add-atom, delete-atom, add-bond, delete-bond, load-structure)
;; so the building functions don't need a lot of optimization.

;; Functions with names like '???-force' run during compute-forces, and
;; they compute the force contributions for each kind of MM2 energy term.
;; These forces are applied to all atoms that participate in computing the
;; energy term. compute-forces runs frequently, so it's worth a lot of
;; effort to optimize the '???-force' functions.

(define (build-length-term m n)
  (entering "build-length-term")
  (let* ((ma (list-ref atom-list m))
	 (na (list-ref atom-list n))
	 (coeffs
	  (lookup-length-coeffs
	   ((ma 'species) 'mm2-index)
	   ((na 'species) 'mm2-index)))
	 (ks (* 0.01 (car coeffs)))
	 (r0 (cadr coeffs)))
    (lambda ()
      (length-force ks r0 ma na m n))))

(define (length-force ks r0 ma na m n)
  (entering "length-force")
  (let* ((rdd (lookup-diff-dist m n))
	 (rd (rdd 'diff))
	 (r (rdd 'distance))
	 (du-dr (* ks (- r r0)))
	 (mult (/ (- du-dr) r)))
    (ma 'add-force (vscale rd mult))
    (na 'add-force (vscale rd (- mult)))))

(define (old-length-force ks r0 ma na)
  (entering "length-force")
  (let* ((rd (vdiff (ma 'position) (na 'position)))
	 (r (safe-sqrt (dot-product rd rd)))
	 (du-dr (* ks (- r r0)))
	 (mult (/ (- du-dr) r)))
    (ma 'add-force (vscale rd mult))
    (na 'add-force (vscale rd (- mult)))))

;; OK, for angles it gets a bit trickier. Let the atom positions be vectors
;; r1, r2, r3, with r2 the vertex. Define the following:
;; L1 = -dotproduct(r1-r2,r3-r2) * (r1-r2) + dotproduct(r1-r2,r1-r2) * (r3-r2)
;; L3 = -dotproduct(r1-r2,r3-r2) * (r3-r2) + dotproduct(r3-r2,r3-r2) * (r1-r2)
;; m1 = L1/|L1|, m3 = L3/|L3|
;; Potential energy is given by U = f(theta), force on atom 1 is
;; f1 = -m1 * dU/dm1 = (-m1/|r1-r2|) * dU/dtheta
;; Likewise f3 = (-m3/|r3-r2|) * dU/dtheta
;; Conservation: f2 = -f1-f3

(define (build-angle-term m n p)
  (entering "build-angle-term")
  (let* ((ma (list-ref atom-list m))
	 (na (list-ref atom-list n))
	 (pa (list-ref atom-list p))
	 (coeffs
	  (lookup-angle-coeffs
	   ((ma 'species) 'mm2-index)
	   ((na 'species) 'mm2-index)
	   ((pa 'species) 'mm2-index)))
	 (kth (car coeffs))
	 (th0 (cadr coeffs)))
    (lambda ()
      (angle-force kth th0 ma na pa m n p))))

(define (angle-force kth th0 ma na pa m n p)
  (entering "angle-force")
  (let* ((dd12 (lookup-diff-dist m n))
	 (dd32 (lookup-diff-dist p n))
	 (th (v-angle dd32 dd12))
	 (tdif (- th th0))
	 (du-dth (* kth (* tdif (+ 1 (* 1.508 tdif tdif)))))
	 (r1r2 (dd12 'diff))
	 (r3r2 (dd32 'diff))
	 (L1 (vdiff (vscale r3r2 (dot-product r1r2 r1r2))
		    (vscale r1r2 (dot-product r1r2 r3r2))))
	 (f1 (vscale L1 (/ du-dth (* (vlen L1) (vlen r1r2)))))
	 (L3 (vdiff (vscale r1r2 (dot-product r3r2 r3r2))
		    (vscale r3r2 (dot-product r1r2 r3r2))))
	 (f3 (vscale L3 (/ du-dth (* (vlen L3) (vlen r3r2)))))
	 (f2 (vscale (vplus f1 f3) -1.0)))
    (ma 'add-force f1)
    (na 'add-force f2)
    (pa 'add-force f3)))

;; To think about torsions, think of the projection of the whole thing into
;; the plane perpendicular to the torqued bond.

(define torsion-whine #f)
(define use-torsion-forces #f)

(define (build-torsion-term m n p q)
  (entering "build-torsion-term")
  (let* ((ma (list-ref atom-list m))
	 (na (list-ref atom-list n))
	 (pa (list-ref atom-list p))
	 (qa (list-ref atom-list q))
	 (coeffs
	  (lookup-torsion-coeffs
	   ((ma 'species) 'mm2-index)
	   ((na 'species) 'mm2-index)
	   ((pa 'species) 'mm2-index)
	   ((qa 'species) 'mm2-index)))
	 (v1 (car coeffs))
	 (v2 (cadr coeffs))
	 (v3 (caddr coeffs)))
    (lambda ()
      (torsion-force v1 v2 v3 ma na pa qa m n p q))))

(define (torsion-force v1 v2 v3 ma na pa qa m n p q)
  (entering "torsion-force")
  (if use-torsion-forces
      (set! torsion-whine #t)))

;; Actually, torsions appear to contribute a negligible amount of force
;; compared even to van der Waals, but they require this atrociously complex
;; calculation. I'm thinking about just not bothering to compute them at all.
;; I think the math here is correct, just woefully inefficient.

;; torsion is broken for the time being
(define (do-not-use-torsion-force v1 v2 v3 ma na pa qa m n p q)
  (let* ((ddmn (lookup-diff-dist m n))
	 (ddpn (lookup-diff-dist p n))
	 (ddqp (lookup-diff-dist q p))
	 (pv (ddmn 'diff))
	 (qv (ddpn 'diff))
	 (rv (ddqp 'diff))
	 (pp (ddmn 'distance))
	 (qq (ddpn 'distance))
	 (rr (ddqp 'distance))
	 (pq (dot-product pv qv))
	 (qr (dot-product qv rv))
	 (alpha (safe-sqrt (/ (* pq pq) qq)))
	 (beta (safe-sqrt (* qq qq)))
	 (gamma (safe-sqrt (/ (* qr qr) qq)))
	 (vm1 (cross-product qv pv))
	 (vq1 (cross-product rv qv))
	 (vm2 (vscale vm1 (/ 1.0 (* (vlen vm1)))))
	 (vq2 (vscale vq1 (/ 1.0 (* (vlen vq1)))))
	 (w (safe-acos (dot-product vm2 vq2)))
	 (du-dw (* -0.0005 (+ (* v1 (sin w))
			      (* -2 v2 (sin (* 2 w)))
			      (* 3 v3 (sin (* 3 w))))))
	 (fm (vscale vm2 (/ du-dw
			    (safe-sqrt (- pp (/ (* pq pq) qq))))))
	 (fq (vscale vq2 (/ du-dw
			    (safe-sqrt (- rr (/ (* qr qr) qq)))))))
    (ma 'add-force fm)
    (qa 'add-force fq)
    (pa 'add-force
	(vdiff (vscale fm (/ alpha beta))
	       (vscale fq (/ (+ gamma beta) beta))))
    (na 'add-force
	(vdiff (vscale fq (/ gamma beta))
	       (vscale fm (/ (+ alpha beta) beta))))))

;; vdw is similar to length force
;; du/dr = 6*evdw*[(r/rvdw)^-7 - (r/rvdw)^-13]
;; Don't forget the factor of 0.001

(define use-vdw-forces #t)

(define (build-vdw-term m n)
  (entering "build-vdw-term")
  (let* ((ma (list-ref atom-list m))  ;; members of atom-list
	 (na (list-ref atom-list n))
	 (evdw (* 0.5 (+ ((ma 'species) 'evdw)
			 ((na 'species) 'evdw))))
	 (rvdw (+ ((ma 'element) 'rvdw)
		  ((na 'element) 'rvdw))))
    (lambda ()
      (vdw-force evdw rvdw ma na m n))))

(define (vdw-force evdw rvdw ma na m n)
  (entering "vdw-force")
  (if use-vdw-forces
      (let* ((rdd (lookup-diff-dist m n))
	     (rd (rdd 'diff))
	     (r (rdd 'distance))
	     (rsq_recip (/ (* rvdw rvdw) (* r r))) ;; (r/rvdw)^-2
	     (r_recip (safe-sqrt rsq_recip))
	     (r6 (* rsq_recip rsq_recip rsq_recip))    ;; (r/rvdw)^-6
	     (du-dr (* 0.006 evdw r_recip r6 (- 1 (* 2 r6))))
	     (mult (/ (- du-dr) r)))
	(ma 'add-force (vscale rd mult))
	(na 'add-force (vscale rd (- mult))))))

;; We only need to run setup-terms when the structure changes: add-atom,
;; delete-atom, add-bond, delete-bond, or load-structure. This flag keeps
;; track of all that.
(define need-to-resetup-terms #t)

(define (compute-forces)
  (entering "compute-forces")
  (if need-to-resetup-terms
      (let ()
	(setup-terms)
	(set! need-to-resetup-terms #f)))
  ;; precalculate difference vectors and distances for all pairs of atoms
  (build-diff-dist-table)
  ;; set all forces to zero, one force vector for each atom
  (dolist (a atom-list) (a 'zero-force))
  ;; if the user defined external forces, put them in now
  (if external-forces (external-forces))
  ;; compute force contributions for each energy term
  (dolist (x term-list) (x)))

;; ******************************************************************
;; Building up a term list from an atom list and bond list

(define (other-end bond atm)
  (cond ((= atm (bond 'first)) (bond 'second))
	((= atm (bond 'second)) (bond 'first))
	(else #f)))

(define whine-about-bond-count #f)

(define (setup-terms)
  (entering "setup-terms")
  ;; for each atom, count its bonds, verify against valence, and reassign
  ;; its mm2 table index according to hybridization
  (let ((n 0) (bc #f))
    (dolist (atom atom-list)
	    (set! bc (make-bond-count n))
	    (if (not (= ((atom 'element) 'total-bonds)
			(bc 'total-bonds)))
		(set! whine-about-bond-count #t))
	    (hybridize-atom atom bc)
	    (set! n (+ n 1))))

  (set! term-list '())

  ;; enumerate length terms
  (dolist (x bond-list)
	  (set! term-list
		(cons (build-length-term (x 'first) (x 'second))
		      term-list)))

  ;; enumerate angle terms
  (do ((AL atom-list (cdr AL))
       (x 0 (+ x 1)))
      ((null? AL))
    (dolist
     (BL bond-list)
     (let ((y (other-end BL x)))
       (if y
	   (dolist (B2L bond-list)
		   (let ((z (other-end B2L y)))
		     (if (and z (> z x))
			 (set! term-list
			       (cons (build-angle-term x y z)
				     term-list)))))))))

  ;; enumerate the torsion terms
  (do ((AL atom-list (cdr AL))
       (w 0 (+ w 1)))
      ((null? AL))
    (do ((BL bond-list (cdr BL)))
	((null? BL))
      (let ((x (other-end (car BL) w)))
	(if x
	    (do ((B2L bond-list (cdr B2L)))
		((null? B2L))
	      (let ((y (other-end (car B2L) x)))
		(if (and y (not (= w y)))
		    (do ((B3L bond-list (cdr B3L)))
			((null? B3L))
		      (let ((z (other-end (car B3L) y)))
			(if (and z (not (= z x)) (> z w))
			    (set! term-list
				  (cons (build-torsion-term w x y z)
					term-list))))))))))))

  ;; enumerate the van der Waals terms (unbonded atom pairs)
  (do ((AL atom-list (cdr AL))
       (x 0 (+ x 1)))
      ((null? AL))
    (do ((A2L (cdr AL) (cdr A2L))
	 (y (+ x 1) (+ y 1))
	 (flag #t))
	((null? A2L))
      (set! flag #t)
      (do ((BL bond-list (cdr BL)))
	  ((null? BL))
	(let ((p (other-end (car BL) x)))
	  (if (and p (= p y))
	      (let ()
		(set! flag #f)
		(set! BL '(4))))))
      ;; Why '(4), you ask? So that after we take the cdr, we get
      ;; BL being '(), the proper terminating condition.
      (if flag
	  (set! term-list
		(cons (build-vdw-term x y)
		      term-list))))))

;; ******************************************************************
;; Potential energy minimization

;; We compute force vectors, and then scale them until to be displacement
;; vectors, where the largest displacement has a magnitude equal to
;; emin-factor (which is in angstroms).

;; currently, these represent angstroms
(define coarse-emin-factor 0.2)
(define fine-emin-factor 0.01)
(define negligible-force 0.0001)

(define emin-factor fine-emin-factor)

(define (emin-step)
  (entering "emin-step")
  (force-rotate-mode)
  (set! whine-about-bond-count #f)
  (set! torsion-whine #f)
  (dotimes
   (i 20)
   (compute-forces)
   (dolist (a atom-list)
	   (let ((f (a 'force)))
	     (if (> (vlen f) negligible-force)
		 (a 'add-pos
		    (vscale f (/ emin-factor (vlen f))))))))
  (if whine-about-bond-count
      (error-msg "Mismatch between bonds and valences"))
  (if torsion-whine
      (warning-msg "Torsion forces temporarily disabled")))
Modified: Thu Jan 23 17:00:00 1997 GMT
Page accessed 6259 times since Sat Apr 17 22:02:57 1999 GMT