|
wolffd@0
|
1 function pot = extension_pot(oldpot, ddom_u, dsizes, ctaildom_u, csizes)
|
|
wolffd@0
|
2 % EXTENSION_POT Extense a stable CG potential.
|
|
wolffd@0
|
3 % pot = extension_pot(oldpot, ddom_u, ctaildom_u, dsizes, csizes)
|
|
wolffd@0
|
4 % ddom_u Added discrete nodes
|
|
wolffd@0
|
5 % ctaildom_u Added continuous tail nodes
|
|
wolffd@0
|
6 % csizes is the size of the tail nodes.
|
|
wolffd@0
|
7
|
|
wolffd@0
|
8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
wolffd@0
|
9 % A CG potential can be extended by adding discrete variables to its %
|
|
wolffd@0
|
10 % domain of continuous variables to its tail %
|
|
wolffd@0
|
11 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
wolffd@0
|
12 ddom = myunion(oldpot.ddom, ddom_u);
|
|
wolffd@0
|
13 ctaildom = myunion(oldpot.ctaildom, ctaildom_u);
|
|
wolffd@0
|
14 cheaddom = oldpot.cheaddom;
|
|
wolffd@0
|
15 udom = myunion(ddom_u, ctaildom_u);
|
|
wolffd@0
|
16 domain = myunion(oldpot.domain, udom);
|
|
wolffd@0
|
17
|
|
wolffd@0
|
18 ns = zeros(1,max(domain));
|
|
wolffd@0
|
19 ns(ddom_u) = dsizes;
|
|
wolffd@0
|
20 ns(ctaildom_u) = csizes;
|
|
wolffd@0
|
21 ns(oldpot.ddom) = oldpot.dsizes;
|
|
wolffd@0
|
22 ns(oldpot.cheaddom) = oldpot.cheadsizes;
|
|
wolffd@0
|
23 ns(oldpot.ctaildom) = oldpot.ctailsizes;
|
|
wolffd@0
|
24
|
|
wolffd@0
|
25 dsizes = ns(ddom);
|
|
wolffd@0
|
26 dsize = prod(ns(ddom));
|
|
wolffd@0
|
27 cheadsizes = ns(cheaddom);
|
|
wolffd@0
|
28 cheadsize = sum(ns(cheaddom));
|
|
wolffd@0
|
29 ctailsizes = ns(ctaildom);
|
|
wolffd@0
|
30 ctailsize = sum(ns(ctaildom));
|
|
wolffd@0
|
31
|
|
wolffd@0
|
32 BZ = zeros(cheadsize, ctailsize);
|
|
wolffd@0
|
33 potarray = cell(1, dsize);
|
|
wolffd@0
|
34 mask = find_equiv_posns(oldpot.ddom, ddom);
|
|
wolffd@0
|
35
|
|
wolffd@0
|
36 tmask = find_equiv_posns(oldpot.ctaildom, ctaildom);
|
|
wolffd@0
|
37 tu = block(tmask, ctailsizes);
|
|
wolffd@0
|
38
|
|
wolffd@0
|
39 for i=1:dsize
|
|
wolffd@0
|
40 sub1 = ind2subv(dsizes, i);
|
|
wolffd@0
|
41 sub2 = sub1(mask);
|
|
wolffd@0
|
42 ind = subv2ind(oldpot.dsizes, sub2);
|
|
wolffd@0
|
43 if isempty(ind)
|
|
wolffd@0
|
44 ind = 1;
|
|
wolffd@0
|
45 end
|
|
wolffd@0
|
46 potc = struct(oldpot.scgpotc{ind});
|
|
wolffd@0
|
47 p = potc.p;
|
|
wolffd@0
|
48 B = BZ;
|
|
wolffd@0
|
49 if ~isempty(B)
|
|
wolffd@0
|
50 B(:, tu) = potc.B;
|
|
wolffd@0
|
51 end
|
|
wolffd@0
|
52 potarray{i} = scgcpot(cheadsize, ctailsize, p, potc.A, B, potc.C);
|
|
wolffd@0
|
53 end
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 pot = scgpot(ddom, cheaddom, ctaildom, ns,potarray);
|
