CNS data reading script

;;;; Copyright 2006 Joel Bard
;;;; Copyright 2006 Paul Emsley
    
  
;; Read in cns coeff-data (given filenames) and a pdb molecule filename to make
;; maps.
;; 
(define (cns->coot 2fofc-coeffs fofc-coeffs model-pdb)


  ;; remove any trailing spaces of string
  (define (trim-trailing-spaces s)
    
    (let f ((ls (reverse (string->list s))))
      
      (cond 
       ((null? ls) "")
       ((eq? #\space (car ls)) (f (cdr ls)))
       (else
	(list->string (reverse ls))))))


  ;; Return '() on bad, else return a 6 membered list.
  ;; 
  (define (get-cell-from-pdb pdb-file)
    
    (let ((s (run-command/strings "awk" (list "/CRYST1/ {printf \"%.3f %.3f %.3f %.3f %.3f %.3f\\n\", $2, $3, $4, $5, $6, $7}" pdb-file) '())))
      
      (if (not (= (length s) 1))
	  '()
	  (car s))))

  ;; return "" on bad, else return a string that is the space group.
  ;; 
  ;; #have to go through all sorts of hoops because sometimes we might
  ;; #get P 21 and sometimes we might get P 1 21 1 for example
  ;; #sftools needs P21, can't handle P1211
  ;; #that being said - there must be a better way
  ;; set SYMM1=`awk 'BEGIN { FIELDWIDTHS="55 11 4" } /CRYST1/ {gsub(" ","",$2); printf "%s\n", $2}' $3`
  ;; 
  (define (get-symm-from-pdb-file pdb-file)
    
    (let ((s (run-command/strings "awk" (list "/CRYST1/ {s=substr($0,56,8); gsub(\" \",\"\",s);print s }" pdb-file) '())))

      (if (not (= (length s) 1))
	  ""
	  (car s))))

  (define (get-symm-from-pdb-file-2 pdb-file)
    
    (let ((s (run-command/strings "awk" (list "BEGIN { FIELDWIDTHS=\"55 11 4\" } /CRYST1/ {printf \"%s\\n\", $2}" 
					      pdb-file) '())))
      (if (not (= (length s) 1))
	  ""
	  ;; because this is not shell sript, we have to be careful
	  ;; to remove trailing spaces in the variable.  Csh does
	  ;; this implicitly.
	  (trim-trailing-spaces (car s)))))

  ;; Return a symmetry string.
  ;; on a problem, return a #f.  Use that to exit (elsewhere).
  ;; 
  (define (get-symm-from-symop-lib symm2)

    (let ((e (getenv "CCP4")))
      (if (not (string? e))
	  (begin 
	    (format #t "CCP4 not set up.~%")
	    #f)
	  (let ((s (run-command/strings "awk" (list (string-append "/" symm2 "/ { print $4}")
						    (append-dir-file 
						     (append-dir-dir (append-dir-dir e "lib") "data")
						     "symop.lib")) 
					'())))
	    (if (not (= (length s) 1))
		#f
		(car s))))))

    ;; main body
    (let* ((map1-prefix (strip-extension 2fofc-coeffs))
	   (map2-prefix (strip-extension  fofc-coeffs))
	   (map1-tmp (string-append map1-prefix "_tmp.pdb"))
	   (map2-tmp (string-append map2-prefix "_tmp.pdb"))
	   (cell (get-cell-from-pdb model-pdb)))
      
      (format #t "map1-prefix: ~s~%" map1-prefix)
      (format #t "map2-prefix: ~s~%" map2-prefix)
      (format #t "cell: ~s~%" cell)
      
      (let ((symm1 (get-symm-from-pdb-file model-pdb))
	    (pdbset-log "cns2coot-pdbset-tmp.log"))
	
	(format #t "symm1: ~s~%" symm1)
	(goosh-command "pdbset" (list "XYZIN" model-pdb "XYZOUT" map1-tmp) 
		       (list (string-append "CELL " cell)
			     (string-append "SPACEGROUP " symm1))
		       pdbset-log #t)
	(let* ((symm2 (get-symm-from-pdb-file-2 map1-tmp))
	       (nov (format #t "symm2: ~s~%" symm2))
	       (symm (get-symm-from-symop-lib symm2)))
	  
	  (if (eq? #f symm)
	      (format #t "Failed to find symm in symop.lib~%")
	      (begin 
		(format #t "INFO:: SYMM is ~s~%" symm)
		(let ((map1-mtz (string-append map1-prefix ".mtz"))
		      (map2-mtz (string-append map2-prefix ".mtz"))
		      (map-coot-mtz (string-append map1-prefix "-coot.mtz"))
		      (cad-log (string-append "cns2coot-cad-tmp.log"))
		      (sftools-1-log  (string-append "cns2coot-sftools-1-tmp.log"))
		      (sftools-2-log  (string-append "cns2coot-sftools-2-tmp.log")))
		  
		  (if (file-exists? map1-mtz)
		      (delete-file map1-mtz))
		  (if (file-exists? map2-mtz)
		      (delete-file map2-mtz))
		  
		  (goosh-command "sftools" '() (list (string-append "read " 2fofc-coeffs)
						     "cns"
						     cell
						     symm
						     "END"
						     "W"
						     "P"
						     "R"
						     "SET LABELS"
						     "FOM"
						     "PHIC"
						     "SCALE"
						     "FWT"
						     "PHWT"
						     (string-append "WRITE " map1-mtz))
				 sftools-1-log #t)
		  
		  (goosh-command "sftools" '() (list (string-append "read " fofc-coeffs)
						     "cns"
						     cell
						     symm
						     "END"
						     "W"
						     "P"
						     "R"
						     "SET LABELS"
						     "FOM"
						     "PHIC"
						     "SCALE"
						     "DELFWT"
						     "PHDELWT"
						     (string-append "WRITE " map2-mtz))
				 sftools-2-log #t)
		  
		  (goosh-command "cad" (list "HKLIN1" map1-mtz 
					     "HKLIN2" map2-mtz 
					     "HKLOUT" map-coot-mtz)
				 (list "LABIN FILE_NUMBER 1 E1=FOM E2=PHIC E3=FWT E4=PHWT"
				       "LABIN FILE_NUMBER 2 E1=DELFWT E2=PHDELWT"
				       "END") cad-log #t)
		  
		  (if (file-exists? map1-mtz)
		      (delete-file map1-mtz))
		  (if (file-exists? map2-mtz)
		      (delete-file map2-mtz))

		  ;; now load them into Coot
		  ;; 
		  (read-pdb model-pdb)
		  (make-and-draw-map-with-reso-with-refmac-params map-coot-mtz "FWT" "PHWT" "" 0 0 0 
								  "Fobs:None-specified" "SigF:None-specified" 
								  "RFree:None-specified" 0 0 0 -1.00 -1.00)

		  (make-and-draw-map-with-reso-with-refmac-params map-coot-mtz "DELFWT" "PHDELWT" "" 0 1 0
								  "Fobs:None-specified" "SigF:None-specified" 
								  "RFree:None-specified" 0 0 0 -1.00 -1.00)
		  )))))))