DataMuseum.dk

top - metrics - download

    Length: 9984 (0x2700)
    Types: TextFile
    Names: »extxbmatrix«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦58ca399f1⟧ »extbib« 
            └─⟦this⟧ 

TextFile

\f




external
real procedure xbmatrix(Na,BL,S,Rnorm,plan,res);
value BL,Rnorm; integer Na,BL; boolean Rnorm,plan;
integer array S;  zone res;
begin
integer  i, j, k, l, m, n, type, p, q, N3, ibl;
boolean sp, nl;  real tdef;
integer array tail(1:10), A(1:50);  array H(1:1);
zone X((abs Na*3+127)//128*128,1,stderror), B(128*2,2,stderror);
real procedure vector(e,i,j);
   value i,j; integer i,j; array e;
   begin
   real l,r; integer k,Ni,Nj;
   r:= 0; Ni:= (i-1)*3; Nj:= (j-1)*3;
   for k:=1 step 1 until 3 do begin
      e(k):= l:= X(Nj+k)-X(Ni+k); r:= r+ l*l end;
   r:= vector:= 1/sqrt(r);
   for k:= 1 step 1 until 3 do e(k):= e(k)*r
end vector;
real procedure dot(a,b);
   array a,b;
   begin real d; integer k;
   d:= 0; for k:=1 step 1 until 3 do d:= d+ a(k)*b(k);
   dot:= d
end dot;
real procedure cross(a,b,c);
   array a,b,c;
   begin real r; integer k;
   c(1):= a(2)*b(3) - a(3)*b(2);
   c(2):= a(3)*b(1) - a(1)*b(3);
   c(3):= a(1)*b(2) - a(2)*b(1);
   r:= 0; for k:=1 step 1 until 3 do r:= r+ c(k)**2;
   cross:= r
end cross;

comment the cartesian coordinates of the atoms are read in;

N3:= abs Na*3; plan:= true; sp:= false add 32;  nl:= false add 10;
open(X,4,<:xtrack:>,0);  tail(1):= (N3+127)//128*2;
monitor(40,X,0,tail); l:= 0;
outrec(X,N3);
for i:=1 step 1 until Na do begin
for j:=1 step 1 until  3 do begin
   read(in,X(l+j)); repeatchar(in); readchar(in,k);
   if k=10 and j<3 then begin
      for j:=j+1 step 1 until 3 do  X(l+j):= 0 end end;
   plan:= plan and X(l+3) = 0; l:= l+3
end i;
if Na<0 then begin
   array Xa(1:N3);
   Na:= -Na; write(res,<:
STRUCTURAL PARAMETERS::>,nl,1);
   if readstruct(res,Xa,Na) then begin
      tdef:= -1; close(X,true); goto alarm end;
   for i:=1 step 1 until N3 do begin
      X(i):= Xa(i); if i mod 3 = 0 then plan:= plan and X(i)=0
   end;
   write(res,nl,2)
end Na<0;
write(res,<:
CARTESIAN COORDINATES:

Atom no.          x          y:>);
if -,plan then write(res,<:          z:>);
write(res,nl,1);  l:= 0;  n:= if plan then 2 else 3;
for i:=1 step 1 until Na do begin
   write(res,nl,1,<<dddd>,i,<:        :>);
   for j:=1 step 1 until n do 
      write(res,<<-ddd.dddddd>,X(l+j)); l:= l+3
end i;

open(B,4,<:btrack:>,0);  k:= 0;
for ibl:=1 step 1 until BL do
   k:= k + (S(ibl,1)*N3+127)//128;
tail(1):= k; monitor(40,B,0,tail);
close(B,true);
write(res,nl,3,<:
INTERNAL COORDINATES:

type no.    type    subno.        specifications
:>);   open(B,4,<:bitrack:>,0);  tail(1):= N3*2;
monitor(40,B,0,tail);
type:= 0; A(1):= 1;  read(in,p);


comment the types of internal coordinates are read in and the
        B-matrix elements are calculated;

newtype:
readstring(in,H,1); tdef:= H(1);
type:= type+1;  A(type+1):= A(type)+p;

if tdef = real <:comb:> then begin
outrec(B,128);
for j:=1,j+1 while k<>41 do begin
   read(in,B(j));
   repread: repeatchar(in); readchar(in,k);
   if k=47 then begin
      read(in,tdef);  B(j):= B(j)/tdef; goto repread
   end;
   if k=115 then begin
      B(j):= sign(B(j))*sqrt(abs B(j)); goto repread
   end end;
setposition(B,0,A(type+1)-1); A(type):= -A(type); tdef:= 0
end;

if tdef = real <:str:> then begin
array e(1:3);
integer Ni,Nj; real r;
write(res,nl,1,<<dddd>,type,<:         str:>);
for m:=1 step 1 until p do begin
   read(in,i,j);  outrec(B,128);
   write(res,sp,if m=1 then 4 else 20,
         <<dddd>,m,<:    :>,i,j,nl,1);
   for q:=1 step 1 until N3 do B(q):= 0;
   Ni:= (i-1)*3;  Nj:= (j-1)*3;   r:= vector(e,i,j);
   for q:=1 step 1 until 3 do begin
      if Rnorm then e(q):= e(q)*r;
      B(Ni+q):= -e(q);  B(Nj+q):= e(q) end
end;  tdef:= 0
end;

if tdef = real <:bend:> then begin
integer Ni,Nj,Nk;
real cos,ri,rk,a,b;
array ei,ek(1:3);
write(res,nl,1,<<dddd>,type,<:        bend:>);
for m:=1 step 1 until p do begin
   read(in,i,j,k);  outrec(B,128);
   write(res,sp,if m=1 then 4 else 20,
         <<dddd>,m,<:    :>,i,j,k,nl,1);
   for q:=1 step 1 until N3 do B(q):= 0;
   Ni:= (i-1)*3; Nj:= (j-1)*3; Nk:= (k-1)*3;
   ri:= vector(ei,j,i); rk:= vector(ek,j,k); cos:= dot(ei,ek);
   a:= 1/sqrt(1-cos*cos);
   ri:= ri*a;     rk:= rk*a;
   for q:=1 step 1 until 3 do begin
      B(Ni+q):= a:= (ei(q)*cos- ek(q))*ri;
      B(Nk+q):= b:= (ek(q)*cos- ei(q))*rk;
      B(Nj+q):= -a-b end
end;  tdef:= 0
end;

if tdef = real <:outo:> then begin
integer Ni,Nj,Nk,Nl;
real ri,rj,rl,sin,csc,cot,a,b,c;
array ei,ej,el,cij,cjl,cli(1:3);
write(res,nl,1,<<dddd>,type,<:        outo:>);
for m:=1 step 1 until p do begin
   read(in,i,j,k,l);  outrec(B,128);
   write(res,sp,if m=1 then 4 else 20,
         <<dddd>,m,<:    :>,i,j,k,l,nl,1);
   for q:=1 step 1 until N3 do B(q):= 0;
   Ni:= (i-1)*3; Nj:= (j-1)*3; Nk:= (k-1)*3; Nl:= (l-1)*3;
   ri:= vector(ei,k,i); rj:= vector(ej,k,j); rl:= vector(el,k,l);
   csc:= 1/sqrt(cross(ei,ej,cij));  cross(ej,el,cjl);  cross(el,ei,cli);
   cot:= dot(ei,ej)*csc;  sin:= dot(cij,el)*csc;
   a:= 1/sqrt(1-sin*sin);
   ri:= ri*a*csc;  rj:= rj*a*csc;  rl:= rl*a;
   for q:=1 step 1 until 3 do begin
      B(Ni+q):= a:= (cjl(q)-(ei(q)*csc-ej(q)*cot)*sin)*ri;
      B(Nj+q):= b:= (cli(q)-(ej(q)*csc-ei(q)*cot)*sin)*rj;
      B(Nl+q):= c:= (cij(q)*csc-el(q)*sin)*rl;
      B(Nk+q):= -a-b-c end
end;  tdef:= 0
end;

if tdef = real <:lin:> then begin
integer Ni,Nj,Nk;
real ri,rk,a,b;
array ei,v(1:3);
write(res,nl,1,<<dddd>,type,<:         lin:>);
for m:=1 step 1 until p do begin
   read(in,i,j,k,v);  outrec(B,128);
   write(res,sp,if m=1 then 4 else 20,
         <<dddd>,m,<:    :>,i,j,k);
   for q:= 1 step 1 until 3 do write(res,<<-dd.dd>,v(q));
   write(res,nl,1);
   for q:=1 step 1 until N3 do B(q):= 0;
   Ni:= (i-1)*3; Nj:= (j-1)*3; Nk:= (k-1)*3;
   rk:= vector(ei,j,k); ri:= vector(ei,j,i);
   for q:=1 step 1 until 3 do v(q):= v(q)-X(Nj+q);  a:= dot(v,ei);
   for q:=1 step 1 until 3 do v(q):= v(q)-ei(q)*a;
   a:= 1/sqrt(dot(v,v));  ri:= ri*a;  rk:= rk*a;
   for q:=1 step 1 until 3 do begin
      B(Ni+q):= a:= v(q)*ri;
      B(Nk+q):= b:= v(q)*rk;
      B(Nj+q):= -a-b end
end;  tdef:= 0
end;

if tdef = real <:tors:> then begin
integer Ni,Nj,Nk,Nl;
real r,ri,rl,cscj,csck,cosj,cosk;
array ei,el,ejk,cj,ck(1:3);
write(res,nl,1,<<dddd>,type,<:        tors:>);
for m:=1 step 1 until p do begin
   read(in,i,j,k,l);  outrec(B,128);
   write(res,sp,if m=1 then 4 else 20,
         <<dddd>,m,<:    :>,i,j,k,l,nl,1);
   for q:=1 step 1 until N3 do B(q):= 0;
   Ni:= (i-1)*3; Nj:= (j-1)*3; Nk:= (k-1)*3; Nl:= (l-1)*3;
   ri:= vector(ei,j,i); rl:= vector(el,k,l); r:= vector(ejk,j,k);
   cscj:= 1/cross(ei,ejk,cj);  csck:= 1/cross(ejk,el,ck);
   cosj:= dot(ei,ejk)*r*cscj;  cosk:= -dot(el,ejk)*r*csck;
   cscj:= cscj*ri;             csck:= csck*rl;
   for q:=1 step 1 until 3 do begin
      B(Ni+q):= r:= cj(q)*cscj;   B(Nj+q):= -r+cj(q)*cosj-ck(q)*cosk;
      B(Nl+q):= r:= ck(q)*csck;   B(Nk+q):= -r+ck(q)*cosk-cj(q)*cosj end
end;  tdef:= 0
end;

comment additional types can be included here;

if tdef<>0 then begin
   write(res,<:<10>***type :>,string tdef,<: undefined.:>);
   setposition(B,0,0); monitor(48,B,0,tail); close(B,true);
   goto alarm end;
read(in,p);  if p>0 then goto newtype;
close(X,true);

comment search for combinations starts here;

for q:=1 step 1 until type do begin
p:= A(q);
if p<0 then begin
integer t,p1;
A(q):= p:= abs p; setposition(B,0,p-1); inrec(B,128); t:= B(1);
write(res,nl,1,<<dddd>,q,<:        comb. of type:>,t,nl,1);
p1:= abs A(t+1) - A(t);  p:= abs A(q+1) - p;
begin
array BB(1:p1,1:N3), comb(1:p*p1);
for i:=p*p1 step -1 until 1 do comb(i):= B(i+1);
l:= 0;  setposition(B,0,A(t)-1);
for i:=1 step 1 until p do begin
   write(res,sp,20,<<dddd>,i,sp,4);
   for j:=1 step 1 until p1 do begin
      l:= l+1; write(res,<<-ddd.dd>,comb(l)) end;
   write(res,nl,1) end;
for i:=1 step 1 until p1 do begin
   inrec(B,128);
   for j:=1 step 1 until N3 do BB(i,j):= B(j) end;
setposition(B,0,A(q)-1);  l:= 0;
for k:=1 step 1 until p do begin outrec(B,128);
for j:=1 step 1 until N3 do begin
   B(j):= 0;  m:= l;
   for i:=1 step 1 until p1 do begin
      m:= m+1; B(j):= B(j) + comb(m)*BB(i,j)
   end end j;  l:= l+p1
end k
end end end q;

if Rnorm then write(res,<:
   Stretchings defined relative to bondlength.:>,nl,1);

comment the U-matrix defining the symmetrycoordinates is read in,
and the BS-matrix is calculated and stored on disc;

m:= 0;  write(res,nl,2,<:
SYMMETRY COORDINATES (defined by  S = U * I):

Species    no.    type    coefficients:>,nl,1);
for ibl:=1 step 1 until BL do begin
zone BS((S(ibl,1)*N3+127)//128*128,1,stderror);
write(res,nl,1,<<dddd>,ibl);  k:= 0;  n:= S(ibl,1);
open(BS,4,<:btrack:>,0); setposition(BS,0,m);
outrec(BS,n*N3);  m:= m+(n*N3+127)//128;
for l:=1 step 1 until n do begin
read(in,type); q:= A(type); p:= A(type+1)-q;
write(res,sp,if l=1 then 7 else 11,<<dd>,l,sp,5,type,sp,3);
begin
real f, g;  array U(1:p), BI(1:p,1:N3);
setposition(B,0,q-1);  g:= 0;
for i:=1 step 1 until p do begin
   read(in,f);  U(i):= f;  g:= g + f*f;  inrec(B,128);
   write(res,<<  -d.d>,f);
   for j:=1 step 1 until N3 do BI(i,j):= B(j)
end;
g:= 1/sqrt(g); for i:=1 step 1 until p do U(i):= U(i)*g;
for j:=1 step 1 until N3 do BS(k+j):= sum(U(i)*BI(i,j),i,1,p);
k:= k+N3; write(res,nl,1)
end end l;
close(BS,true)
end ibl;
setposition(B,0,0); monitor(48,B,0,tail);  close(B,true);

alarm: xbmatrix:= tdef;
end;
end
▶EOF◀