atomes 1.1.14
atomes: an atomic scale modeling tool box
Loading...
Searching...
No Matches
dlp_init.c
Go to the documentation of this file.
1/* This file is part of the 'atomes' software
2
3'atomes' is free software: you can redistribute it and/or modify it under the terms
4of the GNU Affero General Public License as published by the Free Software Foundation,
5either version 3 of the License, or (at your option) any later version.
6
7'atomes' is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
8without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
9See the GNU General Public License for more details.
10
11You should have received a copy of the GNU Affero General Public License along with 'atomes'.
12If not, see <https://www.gnu.org/licenses/>
13
14Copyright (C) 2022-2024 by CNRS and University of Strasbourg */
15
22/*
23* This file: 'dlp_init.c'
24*
25* Contains:
26*
27
28 - The functions to initialize the creation of a classical force field
29
30*
31* List of functions:
32
33 int get_position_in_field_atom_from_model_id (int fat, int at);
34 int get_field_atom_id_from_model_id (field_molecule * fmol, int at);
35 int get_fragment_atom_id_from_model_id (field_molecule * fmol, int at);
36 int get_atom_id_in_fragment_from_model_id (int frag, int at);
37 int set_atom_id (field_atom* at, int c, int p, int id);
38 int get_struct_id_from_atom_id (int ids, int * aid);
39 int prepare_field_atom (int i, int j, int k, int l, int m);
40 int test_for_bonds (field_atom* at, field_atom* bt);
41 int prepare_field_struct (int ids, int sid, int yes_no_num, int * aid);
42 int bonds_between_atoms (int n, field_atom* at, field_atom* bt, int a, int b);
43 int test_for_angles (field_atom* at,
44 field_atom* bt,
45 field_atom* ct);
46 int angles_from_bonds (int n,
47 field_atom* at,
48 field_atom* bt,
49 field_atom* ct);
50 int test_for_dihedrals (field_atom* at,
51 field_atom* bt,
52 field_atom* ct,
53 field_atom* dt);
54 int dihedrals_from_angles (int n,
55 field_atom* at,
56 field_atom* bt,
57 field_atom* ct,
58 field_atom* dt);
59 int impropers_inversion (int n, int stru,
60 int at, int bt, int ct, int dt,
61 int a, int b, int c, int d);
62 int coord_sphere_multiplicity (atom * at, int id, gboolean set_atom_id, gboolean in_field_atom);
63 int find_neighbor_loop (int frag, int aid, int fai,
64 field_atom* fat,
65 field_neighbor * ngmb,
66 int sid, gboolean save_it);
67 int setup_atomic_weight (int seq);
68 int init_vdw (gboolean init);
69
70 gboolean was_not_created_struct (int ids, int num, int * aid);
71 gboolean in_bond (int at, int bd[2]);
72 gboolean are_neighbors (field_neighbor * ngb, int at);
73 gboolean are_in_bond (atom ato, int at);
74 gboolean is_numbering_possible (int frag);
75 gboolean is_this_numbering_possible_for_this_atom (int frag, field_neighbor * ngma, int atom);
76 gboolean id_n_fold_atoms_in_fragment (int frag, int limit, int num_ngb, int search_type,
77 int init, field_neighbor * ngma_init);
78 gboolean id_atoms_in_fragment (int frag, int seed);
79
80 gchar * set_field_atom_name (field_atom* ato, field_molecule * mol);
81
82 void init_all_atoms (int i);
83 void init_all_bonds ();
84 void init_all_angles ();
85 void init_all_dihedrals ();
86 void init_all_impropers_inversions (int stru);
87 void get_weight (int seq);
88 void find_atom_id_in_field_molecule ();
89 void clean_field_struct_list (field_struct * stru);
90 void init_all_field_struct (gboolean init);
91 void init_field_molecule (int i);
92 void init_dlpoly_field (classical_field * new_field);
93 void init_lammps_field (classical_field * new_field);
94
95 field_atom* init_field_atom (int id, int type, int at, int nat, int coo, int * list);
96 field_shell * init_field_shell (int id, int ia, int ib);
97 field_constraint * init_field_constraint (int id, int ia, int ib);
98 field_pmf * init_field_pmf (int id, int num[2], int * list[2], float * w[2]);
99 field_rigid * init_field_rigid (int id, int num, int * list);
100 field_tethered * init_field_tethered (int id, int num);
101 field_prop * init_field_prop (int ti, int key, gboolean show, gboolean use);
102 field_struct * init_field_struct (int st, int ai, int an, int * aid);
103 field_nth_body * init_field_nth_body (int bi, int bd, int * na, int ** ma, int ** ba);
104 field_external * init_field_external (int bi);
105 field_neighbor * get_init_neighbor (int a);
106
107 classical_field * create_force_field_data_structure (int ai);
108
109*/
110
111#include "dlp_field.h"
112#include "interface.h"
113
114extern void print_all_field_struct (field_molecule * mol, int str);
115
116int * astr;
117float val;
118extern int * duplicate_int (int num, int * old_val);
119extern float * duplicate_float (int num, float * old_val);
120int id_atom;
121int multi;
122int a_multi;
123
132int get_position_in_field_atom_from_model_id (int fat, int at)
133{
134 int i;
135 field_atom* tmp_ft = get_active_atom (tmp_fmol -> id, fat);
136 for (i=0; i<tmp_ft -> num; i++)
137 {
138 if (tmp_ft -> list[i] == at) return i;
139 }
140 return -1;
141}
142
151int get_field_atom_id_from_model_id (field_molecule * fmol, int at)
152{
153 int i;
154 field_atom* fat = fmol -> first_atom;
155 while (fat)
156 {
157 for (i=0; i< fat -> num; i++)
158 {
159 if (at == fat -> list[i]) return fat -> id;
160 }
161 fat = fat -> next;
162 }
163 return -1;
164}
165
174int get_fragment_atom_id_from_model_id (field_molecule * fmol, int at)
175{
176 int i;
177 field_atom* fat = fmol -> first_atom;
178 while (fat)
179 {
180 for (i=0; i< fat -> num; i++)
181 {
182 if (at == fat -> list[i]) return fat -> list_id[i];
183 }
184 fat = fat -> next;
185 }
186 return -1;
187}
188
197int get_atom_id_in_fragment_from_model_id (int frag, int at)
198{
199 int i, j, k;
200 field_atom* tmp_ft;
201 for (i=0; i< tmp_fmol -> mol -> natoms; i++)
202 {
203 if (tmp_fmol -> atoms_id[i][frag].a > -1 && tmp_fmol -> atoms_id[i][frag].b > -1)
204 {
205 tmp_ft = get_active_atom (tmp_fmol -> id, tmp_fmol -> atoms_id[i][frag].a);
206 j = tmp_fmol -> atoms_id[i][frag].b;
207 k = tmp_ft -> list[j];
208 if (k == at) return i;
209 }
210 }
211 return -1;
212}
213
224int set_atom_id (field_atom* at, int c, int p, int id)
225{
226 if (c == tmp_fmol -> fragments[0])
227 {
228 id_atom ++;
229 tmp_fmol -> atoms_id[id_atom][0].a = id;
230 tmp_fmol -> atoms_id[id_atom][0].b = p;
231 return id_atom;
232 }
233 else
234 {
235 return -1;
236 }
237}
238
247gchar * set_field_atom_name (field_atom* ato, field_molecule * mol)
248{
249 int i, j;
250 gchar * str;
251 field_atom* tmpatf;
252 switch (ato -> type)
253 {
254 case SPEC_ONLY:
255 str = exact_name(tmp_proj -> chemistry -> label[ato -> sp]);
256 break;
257 case OTHER:
258 tmpatf = mol -> first_atom;
259 j = 1;
260 while (tmpatf)
261 {
262 if (tmpatf -> sp == ato -> sp) j ++;
263 tmpatf = tmpatf -> next;
264 }
265 str = exact_name(g_strdup_printf ("%s-%d-%d", tmp_proj -> chemistry -> label[ato -> sp], j, mol -> id+1));
266 break;
267 default:
268 i = ato -> list[0];
269 j = tmp_proj -> atoms[0][i].coord[ato -> type-1];
270 switch (ato -> type)
271 {
272 case TOTAL_COORD_AND_SPEC:
273 str = g_strdup_printf ("%s-%d-%d",
274 exact_name(tmp_proj -> chemistry -> label[ato -> sp]),
275 tmp_coord -> geolist[0][ato -> sp][j],
276 mol -> id+1);
277 break;
278 case PARTIAL_COORD_AND_SPEC:
279 str = g_strdup_printf ("%s-%d-%d",
280 exact_name(tmp_proj -> chemistry -> label[ato -> sp]),
281 j+1,
282 mol -> id+1);
283 /*str = g_strdup_printf ("%s-%d",
284 exact_name (env_name (tmp_proj, j, ato -> sp, ato -> type, NULL)),
285 mol -> id+1);*/
286 break;
287 }
288 }
289 return str;
290}
291
304field_atom* init_field_atom (int id, int type, int at, int nat, int coo, int * list)
305{
306 int h, i, j, k, l, m, n;
307 field_atom* ato;
308 ato = g_malloc0 (sizeof*ato);
309 ato -> id = id;
310 ato -> type = type;
311 ato -> num = nat;
312 ato -> list = allocint (ato -> num);
313 ato -> list_id = allocint (ato -> num);
314 ato -> frozen_id = allocbool (ato -> num);
315 h = 0;
316 switch (type)
317 {
318 case OTHER:
319 for (i=0; i<tmp_fmol -> multi; i++)
320 {
321 for (j=0; j < nat/tmp_fmol -> multi; j++)
322 {
323 k = list[j];
324 l = tmp_fmol -> atoms_id[k][i].a;
325 m = tmp_fmol -> atoms_id[k][i].b;
326 n = get_active_atom (tmp_fmol -> id, l) -> list[m];
327 ato -> list[h] = n;
328 ato -> list_id[h] = k;
329 ato -> frozen_id[h] = FALSE;
330 tmp_fmol -> atoms_id[k][i].a = ato -> id;
331 tmp_fmol -> atoms_id[k][i].b = h;
332 tmp_proj -> atoms[0][n].faid = ato -> id;
333 h ++;
334 }
335 }
336 break;
337 default:
338 for (i=0; i < tmp_fmol -> multi; i++)
339 {
340 j = tmp_fmol -> fragments[i];
341 for (k=0; k < tmp_proj -> natomes; k++)
342 {
343 if (tmp_proj -> atoms[0][k].coord[2] == j && tmp_proj -> atoms[0][k].sp == at)
344 {
345 if (type == SPEC_ONLY)
346 {
347 ato -> list[h] = k;
348 ato -> frozen_id[h] = FALSE;
349 ato -> list_id[h] = set_atom_id (ato, tmp_proj -> atoms[0][k].coord[2], h, id);
350 h ++;
351 }
352 else
353 {
354 if (tmp_proj -> atoms[0][k].coord[type-1] == coo)
355 {
356 ato -> list[h] = k;
357 ato -> frozen_id[h] = FALSE;
358 ato -> list_id[h] = set_atom_id (ato, tmp_proj -> atoms[0][k].coord[2], h, id);
359 h ++;
360 }
361 }
362 }
363 }
364 }
365 break;
366 }
367#ifdef DEBUG
368 if (h != ato -> num) g_debug ("Something is wrong in with this field atom: id= %d, sp= %d, co= %d", id, at, coo);
369#endif
370 ato -> sp = at;
371
372 ato -> mass = 0.0;//get_force_field_atom_mass (at, 0);
373 ato -> fid = (ato -> mass == 0.0) ? -1 : tmp_field -> type;
374 ato -> afid = (ato -> mass == 0.0) ? -1 : 0;
375 ato -> mass = (ato -> mass == 0.0) ? tmp_proj -> chemistry -> chem_prop[CHEM_M][at] : ato -> mass;
376 ato -> name = g_strdup_printf ("%s", set_field_atom_name (ato, tmp_fmol));
377 ato -> charge = 0.0;
378 ato -> frozen = 0;
379 ato -> show = FALSE;
380 ato -> next = NULL;
381 ato -> prev = NULL;
382 return ato;
383}
384
394field_shell * init_field_shell (int id, int ia, int ib)
395{
396 field_shell * shell;
397 shell = g_malloc0 (sizeof*shell);
398 shell -> id = id;
399 shell -> ia[0] = ia;
400 shell -> ia[1] = ib;
401 shell -> k2 = 0.0;
402 shell -> k4 = 0.0;
403 shell -> show = FALSE;
404 shell -> use = TRUE;
405 shell -> next = NULL;
406 shell -> prev = NULL;
407 return shell;
408}
409
419field_constraint * init_field_constraint (int id, int ia, int ib)
420{
421 field_constraint * cons;
422 cons = g_malloc0 (sizeof*cons);
423 cons -> id = id;
424 cons -> ia[0] = ia;
425 cons -> ia[1] = ib;
426 cons -> av = -1.0;
427 cons -> length = 0.0;
428 cons -> show = FALSE;
429 cons -> use = TRUE;
430 cons -> next = NULL;
431 cons -> prev = NULL;
432 return cons;
433}
434
445field_pmf * init_field_pmf (int id, int num[2], int * list[2], float * w[2])
446{
447 field_pmf * pmf;
448 pmf = g_malloc0 (sizeof*pmf);
449 pmf -> id = id;
450 pmf -> av = -1.0;
451 pmf -> length = 0.0;
452 int i;
453 if (num != NULL)
454 {
455 for (i=0; i<2; i++)
456 {
457 pmf -> num[i] = num[i];
458 pmf -> list[i] = duplicate_int (num[i], list[i]);
459 pmf -> weight[i] = duplicate_float (num[i], w[i]);
460 }
461 }
462 else
463 {
464 for (i=0; i<2; i++)
465 {
466 pmf -> num[i] = 0;
467 pmf -> list[i] = NULL;
468 pmf -> weight[i] = NULL;
469 }
470 }
471 pmf -> show = FALSE;
472 pmf -> use = TRUE;
473 pmf -> next = NULL;
474 pmf -> prev = NULL;
475 return pmf;
476}
477
487field_rigid * init_field_rigid (int id, int num, int * list)
488{
489 field_rigid * rig;
490 rig = g_malloc0 (sizeof*rig);
491 rig -> id = id;
492 rig -> num = num;
493 rig -> list = NULL;
494 if (list != NULL) rig -> list = duplicate_int (num, list);
495 rig -> show = FALSE;
496 rig -> use = TRUE;
497 rig -> next = NULL;
498 rig -> prev = NULL;
499 return rig;
500}
501
510field_tethered * init_field_tethered (int id, int num)
511{
512 field_tethered * tet;
513 tet = g_malloc0 (sizeof*tet);
514 tet -> id = id;
515 tet -> num = num;
516 tet -> key = 0;
517 tet -> val = allocfloat (fvalues[activef][0][0]);
518 tet -> show = FALSE;
519 tet -> use = TRUE;
520 tet -> next = NULL;
521 tet -> prev = NULL;
522 return tet;
523}
524
535field_prop * init_field_prop (int ti, int key, gboolean show, gboolean use)
536{
537 field_prop * prop;
538 prop = g_malloc0 (sizeof*prop);
539 prop -> aid = allocint (struct_id(ti+7));
540 prop -> key = key;
541 prop -> pid = -1;
542 prop -> fpid = -1;
543 prop -> val = allocfloat (fvalues[activef][ti+1][key]);
544 prop -> show = show;
545 prop -> use = use;
546 prop -> next = NULL;
547 prop -> prev = NULL;
548 return prop;
549}
550
559int get_struct_id_from_atom_id (int ids, int * aid)
560{
561 int i, j;
562 field_struct * stru = tmp_fmol -> first_struct[ids];
563 j = struct_id (ids+7);
564 int res = 0;
565 while (stru)
566 {
567 res = stru -> id + 1;
568 for (i=0; i<j; i++)
569 {
570 if (stru -> aid[i] != aid[i])
571 {
572 res = 0;
573 break;
574 }
575 }
576 if (! res)
577 {
578 res = -(stru -> id + 1);
579 for (i=0; i<j; i++)
580 {
581 if (stru -> aid[i] != aid[j-1-i])
582 {
583 res = 0;
584 break;
585 }
586 }
587 }
588 if (res != 0) break;
589 stru = stru -> next;
590 }
591 return res;
592}
593
603gboolean was_not_created_struct (int ids, int num, int * aid)
604{
605 if (num == 0) return TRUE;
606 int j, k;
607 k = struct_id(ids+7);
608 tmp_fstr = tmp_fmol -> first_struct[ids];
609 gboolean res;
610 while (tmp_fstr)
611 {
612 res = FALSE;
613 if (ids < 6)
614 {
615 for (j=0; j<k; j++)
616 {
617 if (tmp_fstr -> aid[j] != aid[j])
618 {
619 res = TRUE;
620 break;
621 }
622 }
623 if (res)
624 {
625 res = FALSE;
626 for (j=0; j<k; j++)
627 {
628 if (tmp_fstr -> aid[k-1-j] != aid[j])
629 {
630 res = TRUE;
631 break;
632 }
633 }
634 }
635 }
636 else if (ids == 6)
637 {
638 for (j=0; j<k; j++)
639 {
640 if (tmp_fstr -> aid[j] != aid[j])
641 {
642 res = TRUE;
643 break;
644 }
645 }
646 if (res)
647 {
648 if (tmp_fstr -> aid[0] != aid[0])
649 {
650 res = TRUE;
651 }
652 else if (tmp_fstr -> aid[3] != aid[3])
653 {
654 res = TRUE;
655 }
656 else if (tmp_fstr -> aid[1] == aid[2] && tmp_fstr -> aid[2] == aid[1])
657 {
658 res = FALSE;
659 }
660 }
661 }
662 else
663 {
664 if (tmp_fstr -> aid[0] != aid[0])
665 {
666 res = TRUE;
667 }
668 else
669 {
670 if (tmp_fstr -> aid[1] != aid[1])
671 {
672 res = TRUE;
673 }
674 else if (tmp_fstr -> aid[2] == aid[2] && tmp_fstr -> aid[3] == aid[3])
675 {
676 res = FALSE;
677 }
678 else if (tmp_fstr -> aid[2] == aid[3] && tmp_fstr -> aid[3] == aid[2])
679 {
680 res = FALSE;
681 }
682 else
683 {
684 res = TRUE;
685 }
686 }
687 }
688 if (! res) return FALSE;
689 tmp_fstr = tmp_fstr -> next;
690 }
691 return TRUE;
692}
693
704field_struct * init_field_struct (int st, int ai, int an, int * aid)
705{
706 field_struct * str;
707 str = g_malloc0 (sizeof*str);
708 str -> st = st;
709 str -> id = ai;
710 if (st == 6 || st == 7)
711 {
712 str -> num = 0;
713 str -> av = 0.0;
714 }
715 else
716 {
717 str -> num = an;
718 str -> av = val;
719 }
720 str -> aid = NULL;
721 if (aid != NULL) str -> aid = duplicate_int (struct_id(st+7), aid);
722 str -> def = init_field_prop (st, 0, FALSE, TRUE);
723 str -> other = NULL;
724 str -> next = NULL;
725 str -> prev = NULL;
726 return str;
727}
728
740field_nth_body * init_field_nth_body (int bi, int bd, int * na, int ** ma, int ** ba)
741{
742 int i, j;
743 field_nth_body * nthbd;
744 nthbd = g_malloc0 (sizeof*nthbd);
745 nthbd -> id = bi;
746 nthbd -> bd = bd;
747 if (! bd)
748 {
749 nthbd -> fpid = allocint(2);
750 nthbd -> fpid[0] = nthbd -> fpid[1] = -1;
751 }
752 j = body_at (bd);
753 // Id in Vdw list
754 if (na == NULL)
755 {
756 nthbd -> na = allocint (j);
757 for (i=0; i<j; i++) nthbd -> na[i] = -1;
758 }
759 else
760 {
761 nthbd -> na = duplicate_int (j, na);
762 }
763 // Molecule id
764 nthbd -> ma = NULL;
765 nthbd -> ma = g_malloc0 (j*sizeof*nthbd -> ma);
766 // Field atom id in molecule
767 nthbd -> a = NULL;
768 nthbd -> a = g_malloc0 (j*sizeof*nthbd -> a);
769 for (i=0; i < j; i++)
770 {
771 if (ma != NULL) nthbd -> ma[i] = duplicate_int (na[i], ma[i]);
772 if (ba != NULL) nthbd -> a[i] = duplicate_int (na[i], ba[i]);
773 }
774 nthbd -> key = 0;
775 nthbd -> val = allocfloat (fvalues[activef][9+bd][nthbd -> key]);
776 nthbd -> show = FALSE;
777 nthbd -> use = TRUE;
778 nthbd -> next = NULL;
779 nthbd -> prev = NULL;
780 return nthbd;
781}
782
790field_external * init_field_external (int bi)
791{
792 field_external * nfext;
793 nfext = g_malloc0 (sizeof*nfext);
794 nfext -> id = bi;
795 nfext -> key = -1;
796 nfext -> val = NULL;
797 nfext -> use = FALSE;
798 nfext -> next = NULL;
799 nfext -> prev = NULL;
800 return nfext;
801}
802
814int prepare_field_atom (int i, int j, int k, int l, int m)
815{
816 if (tmp_fmol -> first_atom == NULL)
817 {
818 tmp_fmol -> first_atom = init_field_atom (i, j, k, l, m, NULL);
819 tmp_fat = tmp_fmol -> first_atom;
820 }
821 else
822 {
823 tmp_fat -> next = init_field_atom (i, j, k, l, m, NULL);
824 tmp_fat -> next -> prev = g_malloc0 (sizeof*tmp_fat -> next -> prev);
825 tmp_fat -> next -> prev = tmp_fat;
826 tmp_fat = tmp_fat -> next;
827 }
828 return 1;
829}
830
838void init_all_atoms (int i)
839{
840 int j, k, l, m, n;
841
842 tmp_fmol -> first_atom = NULL;
843 k = 0;
844 switch (tmp_field -> atom_init)
845 {
846 case SPEC_ONLY:
847 for (j=0; j<tmp_proj -> nspec; j++)
848 {
849 if (tmp_mol -> species[j] > 0)
850 {
851 k += prepare_field_atom (k, tmp_field -> atom_init, j, tmp_mol -> species[j]*tmp_fmol -> multi, 0);
852 }
853 }
854 break;
855 case OTHER:
856
857 break;
858 default:
859 for (j=0; j<tmp_proj -> nspec; j++)
860 {
861 if (tmp_mol -> species[j] > 0)
862 {
863 for (l=0; l < tmp_coord -> totcoord[tmp_field -> atom_init-1]; l++)
864 {
865 n = 0;
866 for (m=0; m < tmp_proj -> natomes; m++)
867 {
868 if (tmp_proj -> atoms[0][m].sp == j
869 && tmp_proj -> atoms[0][m].coord[3] == i
870 && tmp_proj -> atoms[0][m].coord[tmp_field -> atom_init -1] == l) n ++;
871 }
872 if (n > 0)
873 {
874 k += prepare_field_atom (k, tmp_field -> atom_init, j, n, l);
875 }
876 }
877 }
878 }
879 break;
880 }
881 tmp_fmol -> atoms = k;
882}
883
892gboolean in_bond (int at, int bd[2])
893{
894 if (at == bd[0] || at == bd[1])
895 {
896 return TRUE;
897 }
898 else
899 {
900 return FALSE;
901 }
902}
903
912gboolean are_neighbors (field_neighbor * ngb, int at)
913{
914 int i;
915 for (i=0; i<ngb -> num; i++)
916 {
917 if (ngb -> vois[i] == at) return TRUE;
918 }
919 return FALSE;
920}
921
930gboolean are_in_bond (atom ato, int at)
931{
932 int i;
933 for (i=0; i<ato.numv; i++)
934 {
935 if (ato.vois[i] == at) return TRUE;
936 }
937 return FALSE;
938}
939
948int test_for_bonds (field_atom* at, field_atom* bt)
949{
950 int i, j, k, l, m;
951 m = 0;
952 val = 0.0;
953 for (i=0; i < at -> num; i++)
954 {
955 j = at -> list[i];
956 for (k= 0; k < tmp_proj -> atoms[0][j].numv; k++)
957 {
958 l = tmp_proj -> atoms[0][j].vois[k];
959 if (tmp_proj -> atoms[0][l].faid == bt -> id)
960 {
961 m ++;
962 val += distance_3d (& tmp_proj -> cell, 0, & tmp_proj -> atoms[0][j], & tmp_proj -> atoms[0][l]).length;
963 }
964 }
965 }
966 if (m > 0) val /= m;
967 if (at -> id == bt -> id) m /= 2;
968 return m / tmp_fmol -> multi;
969}
970
981int prepare_field_struct (int ids, int sid, int yes_no_num, int * aid)
982{
983 if (yes_no_num > 0)
984 {
985 if (tmp_fmol -> first_struct[ids] == NULL)
986 {
987 tmp_fmol -> first_struct[ids] = init_field_struct (ids, sid, yes_no_num, aid);
988 tmp_fstr = tmp_fmol -> first_struct[ids];
989 return 1;
990 }
991 else if (was_not_created_struct (ids, sid, aid))
992 {
993 tmp_fstr = tmp_fmol -> first_struct[ids];
994 while (tmp_fstr -> next) tmp_fstr = tmp_fstr -> next;
995 tmp_fstr -> next = init_field_struct (ids, sid, yes_no_num, aid);
996 tmp_fstr -> next -> prev = tmp_fstr;
997 tmp_fstr = tmp_fstr -> next;
998 return 1;
999 }
1000 }
1001 return 0;
1002}
1003
1015int bonds_between_atoms (int n, field_atom* at, field_atom* bt, int a, int b)
1016{
1017 if ((at -> sp == a && bt -> sp == b) || (at -> sp == b && bt -> sp == a))
1018 {
1019 astr[0] = at -> id;
1020 astr[1] = bt -> id;
1021 n += prepare_field_struct (0, n, test_for_bonds (at, bt), astr);
1022 }
1023 return n;
1024}
1025
1031void init_all_bonds ()
1032{
1033 int j, k, l;
1034 tmp_fmol -> first_struct[0] = NULL;
1035 k = 0;
1036 for (j=0; j< tmp_proj -> nspec; j++)
1037 {
1038 tmp_fat = tmp_fmol -> first_atom;
1039 while (tmp_fat)
1040 {
1041 tmp_fbt = tmp_fat;
1042 while (tmp_fbt)
1043 {
1044 k = bonds_between_atoms (k, tmp_fat, tmp_fbt, j, j);
1045 tmp_fbt = tmp_fbt -> next;
1046 }
1047 tmp_fat = tmp_fat -> next;
1048 }
1049 }
1050 for (j=0; j<tmp_proj -> nspec-1; j++)
1051 {
1052 for (l=j+1; l<tmp_proj -> nspec; l++)
1053 {
1054 tmp_fat = tmp_fmol -> first_atom;
1055 while (tmp_fat)
1056 {
1057 tmp_fbt = tmp_fat -> next;
1058 while (tmp_fbt)
1059 {
1060 k = bonds_between_atoms (k, tmp_fat, tmp_fbt, j, l);
1061 tmp_fbt = tmp_fbt -> next;
1062 }
1063 tmp_fat = tmp_fat -> next;
1064 }
1065 }
1066 }
1067 tmp_fmol -> nstruct[0] = k;
1068}
1069
1081int test_for_angles (field_atom* at,
1082 field_atom* bt,
1083 field_atom* ct)
1084{
1085 int i, j, k, l, m, n, o;
1086 val = 0.0;
1087 o = 0;
1088 for (i=0; i<at -> num; i++)
1089 {
1090 j = at -> list[i];
1091 for (k=0; k<tmp_proj -> atoms[0][j].numv; k++)
1092 {
1093 l = tmp_proj -> atoms[0][j].vois[k];
1094 if (tmp_proj -> atoms[0][l].numv >= 2 && tmp_proj -> atoms[0][l].faid == bt -> id)
1095 {
1096 for (m=0; m<tmp_proj -> atoms[0][l].numv; m++)
1097 {
1098 n = tmp_proj -> atoms[0][l].vois[m];
1099 if (n != j && tmp_proj -> atoms[0][n].faid == ct -> id)
1100 {
1101 o ++;
1102 val += angle_3d (& tmp_proj -> cell, 0,
1103 & tmp_proj -> atoms[0][j],
1104 & tmp_proj -> atoms[0][l],
1105 & tmp_proj -> atoms[0][n]).angle;
1106 }
1107 }
1108 }
1109 }
1110 }
1111 if (o > 0) val /= o;
1112 if (at -> id == ct -> id) o /= 2;
1113 return o / tmp_fmol -> multi;
1114}
1115
1129int angles_from_bonds (int n, field_atom* at,
1130 field_atom* bt,
1131 field_atom* ct)
1132{
1133 astr[0] = at -> id;
1134 astr[1] = bt -> id;
1135 astr[2] = ct -> id;
1136 n += prepare_field_struct (2, n, test_for_angles (at, bt, ct), astr);
1137 return n;
1138}
1139
1145void init_all_angles ()
1146{
1147 int m, p;
1148
1149 tmp_fmol -> first_struct[2] = NULL;
1150 p = 0;
1151 field_struct * tmp_fst = tmp_fmol -> first_struct[0];
1152 for (m=0; m < tmp_fmol -> nstruct[0]; m++)
1153 {
1154 tmp_fat = get_active_atom (tmp_fmol -> id, tmp_fst -> aid[0]);
1155 tmp_fbt = get_active_atom (tmp_fmol -> id, tmp_fst -> aid[1]);
1156 tmp_fct = tmp_fmol -> first_atom;
1157 while (tmp_fct)
1158 {
1159 p = angles_from_bonds (p, tmp_fat, tmp_fbt, tmp_fct);
1160 p = angles_from_bonds (p, tmp_fbt, tmp_fat, tmp_fct);
1161 tmp_fct = tmp_fct -> next;
1162 }
1163 if (tmp_fst -> next != NULL) tmp_fst = tmp_fst -> next;
1164 }
1165 tmp_fmol -> nstruct[2] = p;
1166}
1167
1181int test_for_dihedrals (field_atom* at,
1182 field_atom* bt,
1183 field_atom* ct,
1184 field_atom* dt)
1185{
1186 int i, j, k, l, m, n, o, p, q;
1187 q = 0;
1188 val = 0.0;
1189 for (i=0; i<at -> num; i++)
1190 {
1191 j = at -> list[i];
1192 for (k=0; k<tmp_proj -> atoms[0][j].numv; k++)
1193 {
1194 l = tmp_proj -> atoms[0][j].vois[k];
1195 if (tmp_proj -> atoms[0][l].faid == bt -> id)
1196 {
1197 for (m=0; m<tmp_proj -> atoms[0][l].numv; m++)
1198 {
1199 n = tmp_proj -> atoms[0][l].vois[m];
1200 if (n != j && tmp_proj -> atoms[0][n].faid == ct -> id)
1201 {
1202 for (o=0; o<tmp_proj -> atoms[0][n].numv; o++)
1203 {
1204 p = tmp_proj -> atoms[0][n].vois[o];
1205 if (p != j && p != l && tmp_proj -> atoms[0][p].faid == dt -> id)
1206 {
1207 q ++;
1208 val += dihedral_3d (& tmp_proj -> cell, 0,
1209 & tmp_proj -> atoms[0][j],
1210 & tmp_proj -> atoms[0][l],
1211 & tmp_proj -> atoms[0][n],
1212 & tmp_proj -> atoms[0][p]).angle;
1213 }
1214 }
1215 }
1216 }
1217 }
1218 }
1219 }
1220 if (q > 0) val /= q;
1221 if (at -> id == dt -> id && bt -> id == ct -> id) q /= 2;
1222 return q / tmp_fmol -> multi;
1223}
1224
1240int dihedrals_from_angles (int n,
1241 field_atom* at,
1242 field_atom* bt,
1243 field_atom* ct,
1244 field_atom* dt)
1245{
1246 astr[0] = at -> id;
1247 astr[1] = bt -> id;
1248 astr[2] = ct -> id;
1249 astr[3] = dt -> id;
1250 n += prepare_field_struct (4, n, test_for_dihedrals (at, bt, ct, dt), astr);
1251 return n;
1252}
1253
1259void init_all_dihedrals ()
1260{
1261 int n, p;
1262
1263 tmp_fmol -> first_struct[4] = NULL;
1264 p = 0;
1265 field_struct * tmp_fst = tmp_fmol -> first_struct[2];
1266 for (n=0; n< tmp_fmol -> nstruct[2]; n++)
1267 {
1268 tmp_fat = get_active_atom (tmp_fmol -> id, tmp_fst -> aid[0]);
1269 tmp_fbt = get_active_atom (tmp_fmol -> id, tmp_fst -> aid[1]);
1270 tmp_fct = get_active_atom (tmp_fmol -> id, tmp_fst -> aid[2]);
1271 tmp_fdt = tmp_fmol -> first_atom;
1272 while (tmp_fdt)
1273 {
1274 p = dihedrals_from_angles (p, tmp_fat, tmp_fbt, tmp_fct, tmp_fdt);
1275 p = dihedrals_from_angles (p, tmp_fct, tmp_fbt, tmp_fat, tmp_fdt);
1276 tmp_fdt = tmp_fdt -> next;
1277 }
1278 if (tmp_fst -> next != NULL) tmp_fst = tmp_fst -> next;
1279 }
1280 tmp_fmol -> nstruct[4] = p;
1281}
1282
1301int impropers_inversion (int n, int stru,
1302 int at, int bt, int ct, int dt,
1303 int a, int b, int c, int d)
1304{
1305 astr[0] = at;
1306 astr[1] = bt;
1307 astr[2] = ct;
1308 astr[3] = dt;
1309 double v;
1310 //g_debug ("Improp:: at= %d, bt= %d, ct= %d, dt= %d", at, bt, ct, dt);
1311 //g_debug ("Improp:: a= %d, b= %d, c= %d, d= %d", a, b, c, d);
1312 n += prepare_field_struct (stru, n, 1, astr);
1313 //g_debug ("tmp_fstr -> id= %d", tmp_fstr -> id);
1314 if (stru == 6)
1315 {
1316 v = dihedral_3d (& tmp_proj -> cell, 0,
1317 & tmp_proj -> atoms[0][b],
1318 & tmp_proj -> atoms[0][c],
1319 & tmp_proj -> atoms[0][a],
1320 & tmp_proj -> atoms[0][d]).angle;
1321 }
1322 else
1323 {
1324 v = inversion_3d (& tmp_proj -> cell, 0,
1325 & tmp_proj -> atoms[0][a],
1326 & tmp_proj -> atoms[0][b],
1327 & tmp_proj -> atoms[0][c],
1328 & tmp_proj -> atoms[0][d]).angle;
1329 }
1330 tmp_fstr -> av += v;
1331 tmp_fstr -> num ++;
1332 //g_debug ("tmp_fstr -> id = %d, tmp_fstr -> num = %d, v= %f, av= %f\n", tmp_fstr -> id, tmp_fstr -> num, v, tmp_fstr -> av);
1333 return n;
1334}
1335
1343void init_all_impropers_inversions (int stru)
1344{
1345 int i, j, k, l, m, n, p;
1346 int * atid, * matid;
1347 atid = allocint (tmp_proj -> coord -> cmax+1);
1348 matid = allocint (tmp_proj -> coord -> cmax+1);
1349 tmp_fmol -> first_struct[stru] = NULL;
1350 p = 0;
1351 for (i=0; i<tmp_fmol -> mol -> natoms; i++)
1352 {
1353 matid[0] = tmp_fmol -> atoms_id[i][0].a;
1354 j = tmp_fmol -> atoms_id[i][0].b;
1355 k = atid[0] = get_active_atom (tmp_fmol -> id, matid[0]) -> list[j];
1356 if ((tmp_proj -> atoms[0][k].numv > 2 && stru == 6) || (tmp_proj -> atoms[0][k].numv == 3 && stru == 7))
1357 {
1358 for (m=0; m<tmp_proj -> atoms[0][k].numv; m++)
1359 {
1360 atid[m+1] = tmp_proj -> atoms[0][k].vois[m];
1361 l = tmp_proj -> atoms[0][k].vois[m];
1362 matid[m+1] = tmp_proj -> atoms[0][l].faid;
1363 }
1364 if (stru == 6)
1365 {
1366 for (l=0; l<tmp_proj -> atoms[0][k].numv-2; l++)
1367 {
1368 for (m=l+1; m<tmp_proj -> atoms[0][k].numv-1; m++)
1369 {
1370 for (n=m+1; n<tmp_proj -> atoms[0][k].numv; n++)
1371 {
1372 p = impropers_inversion (p, stru,
1373 matid[0], matid[l+1], matid[m+1], matid[n+1],
1374 atid[0], atid[l+1], atid[m+1], atid[n+1]);
1375 p = impropers_inversion (p, stru,
1376 matid[0], matid[l+1], matid[n+1], matid[m+1],
1377 atid[0], atid[l+1], atid[n+1], atid[m+1]);
1378 p = impropers_inversion (p, stru,
1379 matid[0], matid[m+1], matid[n+1], matid[l+1],
1380 atid[0], atid[m+1], atid[n+1], atid[l+1]);
1381 }
1382 }
1383 }
1384 }
1385 else
1386 {
1387 p = impropers_inversion (p, stru,
1388 matid[0], matid[1], matid[2], matid[3],
1389 atid[0], atid[1], atid[2], atid[3]);
1390 p = impropers_inversion (p, stru,
1391 matid[0], matid[2], matid[1], matid[3],
1392 atid[0], atid[2], atid[1], atid[3]);
1393 p = impropers_inversion (p, stru,
1394 matid[0], matid[3], matid[1], matid[2],
1395 atid[0], atid[3], atid[1], atid[2]);
1396 }
1397 }
1398 }
1399 tmp_fmol -> nstruct[stru] = p;
1400 tmp_fstr = tmp_fmol -> first_struct[stru];
1401 for (i=0; i<tmp_fmol -> nstruct[stru]; i++)
1402 {
1403 tmp_fstr -> def -> use = FALSE;
1404 tmp_fstr -> av /= tmp_fstr -> num;
1405 if (tmp_fstr -> next != NULL) tmp_fstr = tmp_fstr -> next;
1406 }
1407 g_free (atid);
1408 g_free (matid);
1409}
1410
1420int coord_sphere_multiplicity (atom * at, int id, gboolean set_atom_id)
1421{
1422 int j, k, l, m;
1423 gboolean not_alone;
1424 multi = 0;
1425 a_multi = 0;
1426 tmp_fct = tmp_fmol -> first_atom;
1427 while (tmp_fct)
1428 {
1429 if (tmp_fct -> sp == at -> sp)
1430 {
1431 for (j=0; j<tmp_fct -> num; j++)
1432 {
1433 k = tmp_fct -> list[j];
1434 not_alone = TRUE;
1435 if (k != at -> id)
1436 {
1437 for (l=0; l<2; l++)
1438 {
1439 if (tmp_proj -> atoms[0][k].coord[l] != at -> coord[l])
1440 {
1441 not_alone = FALSE;
1442 break;
1443 }
1444 }
1445 if (not_alone)
1446 {
1447 multi ++;
1448 if (tmp_fct -> id == tmp_fat -> id) a_multi ++;
1449 if (set_atom_id)
1450 {
1451 l = tmp_proj -> atoms[0][k].coord[2];
1452 for (m=1; m<tmp_fmol -> multi; m++)
1453 {
1454 if (l == tmp_fmol -> fragments[m])
1455 {
1456 tmp_fmol -> atoms_id[id][m].a = tmp_fct -> id;
1457 tmp_fmol -> atoms_id[id][m].b = j;
1458 tmp_fct -> list_id[j] = id;
1459 }
1460 }
1461 }
1462 }
1463 }
1464 }
1465 }
1466 tmp_fct = tmp_fct -> next;
1467 }
1468 return multi;
1469}
1470
1471int assigned;
1472int * i_weight; // For the molecule, indice = atom id for the molecule
1473int * n_weight; // For the analyzed multi, indice = atom id for the project
1474field_neighbor * init_ngb = NULL;
1475
1483field_neighbor * get_init_neighbor (int a)
1484{
1485 field_neighbor * ngb = NULL;
1486 ngb = init_ngb;
1487 while (ngb)
1488 {
1489 if (ngb -> id == a) break;
1490 ngb = ngb -> next;
1491 }
1492 return ngb;
1493}
1494
1502gboolean is_numbering_possible (int frag)
1503{
1504 int i, j, k, l, m, n, o , p;
1505 field_neighbor * ngma, * ngmb;
1506 field_atom* fat, * fbt;
1507
1508 for (i=0; i<tmp_fmol -> mol -> natoms-1; i++)
1509 {
1510 if (tmp_fmol -> atoms_id[i][frag].a > -1 && tmp_fmol -> atoms_id[i][frag].b > -1)
1511 {
1512 ngma = get_init_neighbor (i);
1513 j = tmp_fmol -> atoms_id[i][frag].a;
1514 k = tmp_fmol -> atoms_id[i][frag].b;
1515 fat = get_active_atom(tmp_fmol -> id, j);
1516 l = fat -> list[k];
1517 for (m=i+1; m<tmp_fmol -> mol -> natoms; m++)
1518 {
1519 if (tmp_fmol -> atoms_id[m][frag].a > -1 && tmp_fmol -> atoms_id[m][frag].b > -1)
1520 {
1521 ngmb = get_init_neighbor (m);
1522 n = tmp_fmol -> atoms_id[m][frag].a;
1523 o = tmp_fmol -> atoms_id[m][frag].b;
1524 fbt = get_active_atom(tmp_fmol -> id, n);
1525 p = fbt -> list[o];
1526 if (are_in_bond (tmp_proj -> atoms[0][l], p) != are_neighbors (ngma, m))
1527 {
1528 return FALSE;
1529 }
1530 if (are_in_bond (tmp_proj -> atoms[0][p], l) != are_neighbors (ngmb, i))
1531 {
1532 return FALSE;
1533 }
1534 }
1535 }
1536 }
1537 }
1538 return TRUE;
1539}
1540
1550gboolean is_this_numbering_possible_for_this_atom (int frag, field_neighbor * ngma, int atom)
1551{
1552 int m, n, o, p;
1553 field_atom* fbt;
1554
1555 for (m=0; m<tmp_fmol -> mol -> natoms; m++)
1556 {
1557 if (m != ngma -> id)
1558 {
1559 if (tmp_fmol -> atoms_id[m][frag].a > -1 && tmp_fmol -> atoms_id[m][frag].b > -1)
1560 {
1561 n = tmp_fmol -> atoms_id[m][frag].a;
1562 o = tmp_fmol -> atoms_id[m][frag].b;
1563 fbt = get_active_atom(tmp_fmol -> id, n);
1564 p = fbt -> list[o];
1565 if (are_in_bond (tmp_proj -> atoms[0][atom], p) != are_neighbors (ngma, m))
1566 {
1567 return FALSE;
1568 }
1569 }
1570 }
1571 }
1572 return TRUE;
1573}
1574
1591int find_neighbor_loop (int frag, int aid, int fai,
1592 field_atom* fat,
1593 field_neighbor * ngmb,
1594 int sid, gboolean save_it)
1595{
1596 int i, j, k;
1597 int iter = 0;
1598
1599 i = fat -> num / tmp_fmol -> multi;
1600 for (j=frag*i; j < (frag+1)*i; j++)
1601 {
1602 if (fat -> list_id[j] < 0)
1603 {
1604 k = fat -> list[j];
1605 if (are_in_bond (tmp_proj -> atoms[0][k], aid))
1606 {
1607 //g_debug ("Are linked, n_weight[%d]= %d, i_weight[%d]= %d", k+1, n_weight[k], j+1, i_weight[j]);
1608 if (n_weight[k] == i_weight[ngmb -> id])
1609 {
1610 tmp_fmol -> atoms_id[fai][frag].a = fat -> id;
1611 tmp_fmol -> atoms_id[fai][frag].b = j;
1612 fat -> list_id[j] = fai;
1613 if (save_it && iter == sid) return k;
1614 if (is_this_numbering_possible_for_this_atom (frag, ngmb, k)) iter ++;
1615 tmp_fmol -> atoms_id[fai][frag].a = -1;
1616 tmp_fmol -> atoms_id[fai][frag].b = -1;
1617 fat -> list_id[j] = -1;
1618 }
1619 }
1620 }
1621 }
1622 return iter;
1623}
1624
1638gboolean id_n_fold_atoms_in_fragment (int frag, int limit, int num_ngb, int search_type,
1639 int init, field_neighbor * ngma_init)
1640{
1641 int i, j, k, l, m, n, o, p, q, r, s, t;
1642 field_neighbor * ngma, * ngmb, * ngmc;
1643 field_atom* fat, * fbt;
1644 if (assigned < limit)
1645 {
1646 if (search_type)
1647 {
1648 i = ngma_init -> id;
1649 j = tmp_fmol -> atoms_id[i][frag].a;
1650 for (k=0; k<ngma_init -> num; k++)
1651 {
1652 l = ngma_init -> vois[k];
1653 ngma = get_init_neighbor (l);
1654 if (tmp_fmol -> atoms_id[l][frag].a < 0 && tmp_fmol -> atoms_id[l][frag].b < 0)
1655 {
1656 if (ngma -> num > num_ngb)
1657 {
1658 m = tmp_fmol -> atoms_id[l][0].a;
1659 fat = get_active_atom (tmp_fmol -> id, m);
1660
1661 n = find_neighbor_loop (frag, init, l, fat, ngma, 0, FALSE);
1662 for (o=0; o<n; o++)
1663 {
1664 p = find_neighbor_loop (frag, init, l, fat, ngma, o, TRUE);
1665 assigned ++;
1666 if (assigned == limit)
1667 {
1668 return TRUE;
1669 }
1670 else if (id_n_fold_atoms_in_fragment(frag, limit, num_ngb, 1, p, ngma))
1671 {
1672 return TRUE;
1673 }
1674 else
1675 {
1676 for (q=0; q<tmp_fmol -> mol -> natoms; q++)
1677 {
1678 if (q != l && tmp_fmol -> atoms_id[q][frag].a > -1 && tmp_fmol -> atoms_id[q][frag].b > -1)
1679 {
1680 ngmb = get_init_neighbor(q);
1681 t = 0;
1682 for (r=0; r<ngmb -> num; r++)
1683 {
1684 s = ngmb -> vois[r];
1685 ngmc = get_init_neighbor(s);
1686 if (tmp_fmol -> atoms_id[s][frag].a < 0 && tmp_fmol -> atoms_id[s][frag].b < 0 && ngmc -> num > num_ngb) t ++;
1687 }
1688 if (t)
1689 {
1690 r = tmp_fmol -> atoms_id[q][frag].a;
1691 fbt = get_active_atom (tmp_fmol -> id, r);
1692 s = tmp_fmol -> atoms_id[q][frag].b;
1693 t = fbt -> list[s];
1694 if (id_n_fold_atoms_in_fragment(frag, limit, num_ngb, 1, t, ngmb))
1695 {
1696 return TRUE;
1697 }
1698 }
1699 }
1700 }
1701 if (id_n_fold_atoms_in_fragment(frag, limit, num_ngb, 0, -1, NULL))
1702 {
1703 return TRUE;
1704 }
1705 else
1706 {
1707 assigned --;
1708 tmp_fmol -> atoms_id[l][frag].a = -1;
1709 p = tmp_fmol -> atoms_id[l][frag].b;
1710 fat -> list_id[p] = -1;
1711 tmp_fmol -> atoms_id[l][frag].b = -1;
1712 }
1713 }
1714 }
1715 }
1716 }
1717 }
1718 return FALSE;
1719 }
1720 else
1721 {
1722 for (i=0; i<tmp_fmol -> mol -> natoms; i++)
1723 {
1724 if (tmp_fmol -> atoms_id[i][frag].a < 0 && tmp_fmol -> atoms_id[i][frag].b <0)
1725 {
1726 j = 0;
1727 ngma = get_init_neighbor (i);
1728 if (ngma -> num > num_ngb)
1729 {
1730 for (k=0; k<tmp_proj -> natomes; k++)
1731 {
1732 if (tmp_proj -> atoms[0][k].coord[2] == tmp_fmol -> fragments[frag])
1733 {
1734 j ++;
1735 if ((n_weight[k] == i_weight[i]) && (ngma -> num == tmp_proj -> atoms[0][k].numv))
1736 {
1737 if (get_atom_id_in_fragment_from_model_id (frag, k) < 0)
1738 {
1739 assigned ++;
1740 l = get_field_atom_id_from_model_id (tmp_fmol, k);
1741 tmp_fmol -> atoms_id[i][frag].a = l;
1742 m = get_position_in_field_atom_from_model_id (l, k);
1743 tmp_fmol -> atoms_id[i][frag].b = m;
1744 get_active_atom(tmp_fmol -> id, l) -> list_id[m] = i;
1745 if (is_this_numbering_possible_for_this_atom (frag, ngma, k))
1746 {
1747 if (id_n_fold_atoms_in_fragment(frag, limit, num_ngb, 1, k, ngma))
1748 {
1749 return TRUE;
1750 }
1751 else if (id_n_fold_atoms_in_fragment(frag, limit, num_ngb, 0, -1, NULL))
1752 {
1753 return TRUE;
1754 }
1755 }
1756 tmp_fmol -> atoms_id[i][frag].a = -1;
1757 tmp_fmol -> atoms_id[i][frag].b = -1;
1758 get_active_atom(tmp_fmol -> id, l) -> list_id[m] = -1;
1759 assigned --;
1760 }
1761 }
1762 if (j == tmp_fmol -> mol -> natoms)
1763 {
1764 break;
1765 }
1766 }
1767 }
1768 return FALSE;
1769 }
1770 }
1771 }
1772 return FALSE;
1773 }
1774 }
1775 else
1776 {
1777 return TRUE;
1778 }
1779}
1780
1789gboolean id_atoms_in_fragment (int frag, int seed)
1790{
1791 int i, j;
1792 i = 0;
1793 for (j=0; j<tmp_proj -> natomes; j++)
1794 {
1795 if (tmp_proj -> atoms[0][j].coord[2] == tmp_fmol -> fragments[frag])
1796 {
1797 if (tmp_proj -> atoms[0][j].numv > 1) i++;
1798 }
1799 }
1800 if (seed > -1)
1801 {
1802 field_neighbor * ngb = get_init_neighbor (seed);
1803 if (id_n_fold_atoms_in_fragment (frag, i, 1, 0, seed, ngb))
1804 {
1805 return id_n_fold_atoms_in_fragment (frag, tmp_fmol -> mol -> natoms, 0, 0, seed, ngb);
1806 }
1807 else
1808 {
1809 return FALSE;
1810 }
1811 }
1812 else
1813 {
1814 if (id_n_fold_atoms_in_fragment (frag, i, 1, 0, -1, NULL))
1815 {
1816 return id_n_fold_atoms_in_fragment (frag, tmp_fmol -> mol -> natoms, 0, 0, -1, NULL);
1817 }
1818 else
1819 {
1820 return FALSE;
1821 }
1822 }
1823}
1824
1832void get_weight (int seq)
1833{
1834 int h, i, j, k;
1835 for (i=0; i<tmp_proj -> natomes; i++)
1836 {
1837 for (j=0; j<tmp_proj -> atoms[0][i].numv; j++)
1838 {
1839 k = tmp_proj -> atoms[0][i].vois[j];
1840 n_weight[i] += tmp_proj -> atoms[0][k].numv;
1841 }
1842 }
1843 int * weight = duplicate_int (tmp_proj -> natomes, n_weight);
1844 for (h=1; h<seq; h++)
1845 {
1846 for (i=0; i<tmp_proj -> natomes; i++)
1847 {
1848 for (j=0; j<tmp_proj -> atoms[0][i].numv; j++)
1849 {
1850 k = tmp_proj -> atoms[0][i].vois[j];
1851 weight[i] += n_weight[k];
1852 }
1853 }
1854 for (i=0; i<tmp_proj -> natomes; i++)
1855 {
1856 n_weight[i] = weight[i];
1857 }
1858 }
1859 g_free (weight);
1860
1861#ifdef DEBUG
1862 /*for (i=0; i<tmp_proj -> natomes; i++)
1863 {
1864 j = tmp_proj -> atoms[0][i].sp;
1865 g_debug ("N_WEIGHT:: %s i= %d, n_weight[i] = %d", tmp_proj -> chemistry -> label[j], i+1, n_weight[i]);
1866 }*/
1867#endif
1868}
1869
1877int setup_atomic_weight (int seq)
1878{
1879 int h, i, j, k, l;
1880 field_neighbor * ngm;
1881 i_weight = allocint (tmp_fmol -> mol -> natoms);
1882 init_ngb = NULL;
1883 for (i=0; i<tmp_proj -> natomes; i++)
1884 {
1885 if (tmp_proj -> atoms[0][i].coord[2] == tmp_fmol -> fragments[0])
1886 {
1887 if (init_ngb)
1888 {
1889 ngm -> next = g_malloc0(sizeof*ngm);
1890 ngm -> next -> prev = ngm;
1891 ngm = ngm -> next;
1892 }
1893 else
1894 {
1895 init_ngb = g_malloc0(sizeof*init_ngb);
1896 ngm = init_ngb;
1897 }
1898 ngm -> num = tmp_proj -> atoms[0][i].numv;
1899 ngm -> vois = duplicate_int (ngm -> num, tmp_proj -> atoms[0][i].vois);
1900 j = get_atom_id_in_fragment_from_model_id (0, i);
1901 ngm -> id = j;
1902 for (k=0; k<ngm -> num; k++)
1903 {
1904 l = ngm -> vois[k];
1905 ngm -> vois[k] = get_atom_id_in_fragment_from_model_id (0, l);
1906 i_weight[j] += tmp_proj -> atoms[0][l].numv;
1907 }
1908 }
1909 }
1910
1911 int * weight = duplicate_int (tmp_fmol -> mol -> natoms, i_weight);
1912 for (h=1; h<seq; h++)
1913 {
1914 for (i=0; i<tmp_fmol -> mol -> natoms; i++)
1915 {
1916 ngm = get_init_neighbor (i);
1917 if (ngm)
1918 {
1919 for (j=0; j<ngm -> num; j++)
1920 {
1921 k = ngm -> vois[j];
1922 weight[i] += i_weight[k];
1923 }
1924 }
1925 }
1926 for (i=0; i<tmp_fmol -> mol -> natoms; i++)
1927 {
1928 i_weight[i] = weight[i];
1929 }
1930 }
1931
1932#ifdef DEBUG
1933 /*field_atom* fat;
1934 for (i=0; i<tmp_fmol -> mol -> natoms; i++)
1935 {
1936 fat = get_active_atom (tmp_fmol -> id, tmp_fmol -> atoms_id[i][0].a);
1937 j = tmp_fmol -> atoms_id[i][0].b;
1938 k = fat -> list[j];
1939 l = tmp_proj -> atoms[0][k].sp;
1940 g_debug ("I_WEIGHT:: %s i= %d, i_weight[i] = %d", tmp_proj -> chemistry -> label[l], i, i_weight[i]);
1941 }*/
1942#endif
1943
1944 g_free (weight);
1945
1946 int seed = -1;
1947 for (i=0; i<tmp_fmol -> mol -> natoms; i++)
1948 {
1949 seed = i;
1950 for (j=0; j<tmp_fmol -> mol -> natoms; j++)
1951 {
1952 if (j!=i)
1953 {
1954 if (i_weight[i] == i_weight[j])
1955 {
1956 seed = -1;
1957 break;
1958 }
1959 }
1960 }
1961 if (seed > -1) break;
1962 }
1963 get_weight (seq);
1964 return seed;
1965}
1966
1972void find_atom_id_in_field_molecule ()
1973{
1974 int i, j, k, l;
1975 gboolean sym = FALSE;
1976 gboolean done;
1977 int ** coordnum;
1978 int start = 1;
1979 int seed = -1;
1980 assigned = 0;
1981
1982 coordnum = allocdint(tmp_proj -> nspec, tmp_coord -> totcoord[1]);
1983 tmp_fat = tmp_fmol -> first_atom;
1984 while (tmp_fat)
1985 {
1986 if (tmp_fat -> num == tmp_fmol -> multi)
1987 {
1988 i = tmp_fat -> list_id[0];
1989 // Multiplicity = num of atoms ...
1990 // Unic atom, ie. as much as identical molecular fragments
1991 for (j=1; j<tmp_fmol -> multi; j++)
1992 {
1993 for (k=1; k<tmp_fat -> num; k++)
1994 {
1995 l = tmp_fat -> list[k];
1996 if (tmp_proj -> atoms[0][l].coord[2] == tmp_fmol -> fragments[j])
1997 {
1998 tmp_fmol -> atoms_id[i][j].a = tmp_fmol -> atoms_id[i][0].a;
1999 tmp_fmol -> atoms_id[i][j].b = k;
2000 tmp_fat -> list_id[j] = i;
2001 break;
2002 }
2003 }
2004 }
2005 assigned ++;
2006 }
2007 tmp_fat = tmp_fat -> next;
2008 }
2009
2010 for (i=0; i< tmp_fmol -> mol -> natoms; i++)
2011 {
2012 tmp_fat = get_active_atom (tmp_fmol -> id, tmp_fmol -> atoms_id[i][0].a);
2013 done = TRUE;
2014 for (j=0; j<tmp_fat -> num; j++)
2015 {
2016 if (tmp_fat -> list_id[j] < 0)
2017 {
2018 done = FALSE;
2019 break;
2020 }
2021 }
2022 if (! done)
2023 {
2024 j = tmp_fmol -> atoms_id[i][0].b;
2025 k = tmp_fat -> list[j];
2026 l = coord_sphere_multiplicity (& tmp_proj -> atoms[0][k], i, FALSE);
2027 coordnum[tmp_proj -> atoms[0][k].sp][tmp_proj -> atoms[0][k].coord[1]] = l+1;
2028 if (l == tmp_fmol -> multi - 1 || a_multi == tmp_fmol -> multi - 1)
2029 {
2030 coord_sphere_multiplicity (& tmp_proj -> atoms[0][k], i, TRUE);
2031 assigned ++;
2032 }
2033 }
2034 }
2035 for (i=0; i<tmp_proj -> nspec; i++)
2036 {
2037 for (j=0; j<tmp_coord -> totcoord[1]; j++)
2038 {
2039 if (coordnum[i][j] > 1 && coordnum[i][j] % tmp_fmol -> multi != 0) g_debug ("This should not happen");
2040 coordnum[i][j] /= tmp_fmol -> multi;
2041 }
2042 }
2043 for (k=6; k>0; k--)
2044 {
2045 sym = TRUE;
2046 for (i=0; i<tmp_proj -> nspec; i++)
2047 {
2048 for (j=0; j<tmp_coord -> totcoord[1]; j++)
2049 {
2050 if (coordnum[i][j] > 1 && coordnum[i][j] >= k && coordnum[i][j] % k != 0) sym = FALSE;
2051 }
2052 }
2053 if (sym) break;
2054 }
2055 g_free (coordnum);
2056
2057 if (assigned < tmp_fmol -> mol -> natoms)
2058 {
2059#ifdef DEBUG
2060 if ((sym || assigned == 0))
2061 {
2062 g_debug ("Molecule seems to be symmetric, k= %d", k);
2063 g_debug ("Mol -> id: %d", tmp_fmol -> id);
2064 g_debug ("Mol -> multi= %d", tmp_fmol -> multi);
2065 g_debug ("Mol -> field atoms= %d", tmp_fmol -> atoms);
2066 g_debug ("Mol -> atomes= %d", tmp_fmol -> mol -> natoms);
2067 g_debug ("Assigned so far= %d", assigned);
2068 }
2069#endif
2070
2071 // k is a weight factor
2072 // The values of the total coord of each atom's neighbors
2073 // is sum k times over, this is a safety for systems with
2074 // a peculiar symetry, otherwise it could be very very long
2075 k = 3;
2076 j = assigned;
2077 gboolean done = FALSE;
2078 if (n_weight) g_free (n_weight);
2079 n_weight = NULL;
2080 n_weight = allocint (tmp_proj -> natomes);
2081 while (! done && (seed < 0 || (tmp_fmol -> mol -> natoms && k < tmp_fmol -> mol -> natoms)))
2082 {
2083 seed = setup_atomic_weight (k);
2084 if (seed < 0) k ++;
2085 done = TRUE;
2086 }
2087#ifdef DEBUG
2088 g_debug ("Assigning using: seed= %d, w_factor= %d", seed, k);
2089#endif
2090
2091 for (i=start; i<tmp_fmol -> multi; i++)
2092 {
2093#ifdef DEBUG
2094 g_debug (" -> fragment id: %d", i+1);
2095#endif
2096 assigned = j;
2097 while (! id_atoms_in_fragment (i, seed) && k < tmp_fmol -> mol -> natoms)
2098 {
2099#ifdef DEBUG
2100 g_debug ("**********************************");
2101 g_debug ("*** NOT ABLE TO ASSIGN NUMBERS ***");
2102 g_debug ("*** Mol -> Id= %3d ***", tmp_fmol -> id);
2103 g_debug ("*** Mol -> Multi= %3d ***", tmp_fmol -> multi);
2104 g_debug ("*** Fragment= %3d ***", i+1);
2105 g_debug ("*** k= %3d ***", k);
2106 g_debug ("*** assigned= %10d ***", assigned);
2107 g_debug ("**********************************");
2108#endif
2109 k ++;
2110 if (n_weight) g_free (n_weight);
2111 n_weight = NULL;
2112 n_weight = allocint (tmp_proj -> natomes);
2113 seed = setup_atomic_weight (k);
2114 }
2115 for (l=0; l<tmp_proj -> natomes; l++)
2116 {
2117 if (tmp_proj -> atoms[0][l].coord[2] == tmp_fmol -> fragments[i])
2118 {
2119 tmp_proj -> atoms[0][l].faid = get_field_atom_id_from_model_id (tmp_fmol, l);
2120 tmp_proj -> atoms[0][l].fid = get_atom_id_in_fragment_from_model_id (i, l);
2121 }
2122 }
2123#ifdef DEBUG
2124 g_debug (" -> done !");
2125#endif
2126 }
2127 }
2128 if (i_weight) g_free (i_weight);
2129 i_weight = NULL;
2130 if (n_weight) g_free (n_weight);
2131 n_weight = NULL;
2132 if (init_ngb) g_free (init_ngb);
2133 init_ngb = NULL;
2134#ifdef DEBUG
2135 g_debug ("Assigning completed !");
2136#endif
2137}
2138
2146void clean_field_struct_list (field_struct * stru)
2147{
2148 while (stru)
2149 {
2150 if (stru -> next)
2151 {
2152 stru = stru -> next;
2153 g_free (stru -> prev);
2154 }
2155 else
2156 {
2157 g_free (stru);
2158 stru = NULL;
2159 }
2160 }
2161}
2162
2170void init_all_field_struct (gboolean init)
2171{
2172 int i;
2173 for (i=0; i<tmp_proj -> natomes; i++)
2174 {
2175 if (tmp_proj -> atoms[0][i].coord[2] == tmp_fmol -> fragments[0])
2176 {
2177 tmp_proj -> atoms[0][i].faid = get_field_atom_id_from_model_id (tmp_fmol, i);
2178 tmp_proj -> atoms[0][i].fid = get_atom_id_in_fragment_from_model_id (0, i);
2179 }
2180 }
2181 if (tmp_fmol -> multi > 1 && init) find_atom_id_in_field_molecule ();
2182
2183 for (i=0; i<8; i++)
2184 {
2185 tmp_fmol -> nstruct[i] = 0;
2186 if (tmp_fmol -> first_struct[i] && tmp_fmol -> nstruct[i])
2187 {
2188 clean_field_struct_list (tmp_fmol -> first_struct[i]);
2189 }
2190 tmp_fmol -> first_struct[i] = NULL;
2191 }
2192
2193 // Bonds start
2194 astr = g_malloc0 (2*sizeof*astr);
2195 init_all_bonds ();
2196 g_free (astr);
2197
2198 if (tmp_fmol -> nstruct[0] > 0)
2199 {
2200 tmp_fmol -> first_struct[1] = duplicate_field_struct_list (tmp_fmol -> first_struct[0], TRUE);
2201 tmp_fmol -> nstruct[1] = tmp_fmol -> nstruct[0];
2202 }
2203 // Bonds end
2204
2205 // Angle start
2206 astr = g_malloc0 (3*sizeof*astr);
2207 init_all_angles ();
2208 g_free (astr);
2209 if (tmp_fmol -> nstruct[2] > 0)
2210 {
2211 tmp_fmol -> nstruct[3] = tmp_fmol -> nstruct[2];
2212 tmp_fmol -> first_struct[3] = duplicate_field_struct_list (tmp_fmol -> first_struct[2], TRUE);
2213 }
2214 // Angle end
2215
2216 // Dihedral start
2217 astr = g_malloc0 (4*sizeof*astr);
2218 init_all_dihedrals ();
2219
2220 // Torsional restraints
2221 if (tmp_fmol -> nstruct[4] > 0)
2222 {
2223 tmp_fmol -> nstruct[5] = tmp_fmol -> nstruct[4];
2224 tmp_fmol -> first_struct[5] = duplicate_field_struct_list (tmp_fmol -> first_struct[4], TRUE);
2225 }
2226
2227// Improper and inversion angles
2228 for (i=6; i<8; i++) init_all_impropers_inversions (i);
2229 g_free (astr);
2230}
2231
2239void init_field_molecule (int i)
2240{
2241 int j, k;
2242 tmp_fmol -> id = i;
2243 tmp_fmol -> multi = tmp_mol -> multiplicity;
2244 tmp_fmol -> name = g_strdup_printf("MOL-%d", i+1);
2245 tmp_fmol -> fragments = NULL;
2246 tmp_fmol -> fragments = allocint (tmp_fmol -> multi);
2247 for (j=0; j<tmp_fmol -> multi; j++)
2248 {
2249 tmp_fmol -> fragments[j] = tmp_mol -> fragments[j];
2250 }
2251 tmp_fmol -> mol = tmp_mol;
2252
2253 // Atoms start
2254 id_atom = -1;
2255 tmp_fmol -> atoms_id = g_malloc (tmp_mol -> natoms*sizeof*tmp_fmol -> atoms_id);
2256 for (j=0; j<tmp_mol -> natoms; j++)
2257 {
2258 tmp_fmol -> atoms_id[j] = g_malloc0 (tmp_fmol -> multi*sizeof*tmp_fmol -> atoms_id[j]);
2259 for (k=0; k<tmp_fmol -> multi; k++)
2260 {
2261 tmp_fmol -> atoms_id[j][k].a = -1;
2262 tmp_fmol -> atoms_id[j][k].b = -1;
2263 }
2264 }
2265
2266 init_all_atoms (i);
2267 // Atoms end
2268
2269 // setup_field_molecule_neighbors (i, tmp_proj);
2270
2271 tmp_fmol -> shells = 0;
2272 // Field struct
2273 init_all_field_struct (TRUE);
2274
2275 tmp_fmol -> shells = 0;
2276 tmp_fmol -> first_shell = NULL;
2277 tmp_fmol -> constraints = 0;
2278 tmp_fmol -> first_constraint = NULL;
2279 tmp_fmol -> pmfs = 0;
2280 tmp_fmol -> first_pmf = NULL;
2281 tmp_fmol -> rigids = 0;
2282 tmp_fmol -> first_rigid = NULL;
2283 tmp_fmol -> tethered = 0;
2284 tmp_fmol -> first_tethered = NULL;
2285
2286 tmp_fmol -> show = FALSE;
2287 tmp_fmol -> show_id = FALSE;
2288 tmp_fmol -> next = NULL;
2289 tmp_fmol -> prev = NULL;
2290}
2291
2299int init_vdw (gboolean init)
2300{
2301 int nvdw;
2302 int i, j, k, l;
2303 tmp_fmol = tmp_field -> first_molecule;
2304 i = 0;
2305 while (tmp_fmol)
2306 {
2307 i += tmp_fmol -> atoms;
2308 tmp_fmol = tmp_fmol -> next;
2309 }
2310 nvdw = i;
2311 gchar ** to_be_vdw = g_malloc (nvdw*sizeof*to_be_vdw);
2312 int * vdw_mlist = allocint (nvdw);
2313 int ** vdw_aids = allocdint (nvdw,nvdw);
2314 int ** vdw_mids = allocdint (nvdw,nvdw);
2315 int vdw_na[2];
2316 int * vdw_a[2], * vdw_ma[2];
2317 gboolean add_vdw;
2318 i = 0;
2319 tmp_fmol = tmp_field -> first_molecule;
2320 while (tmp_fmol)
2321 {
2322 tmp_fat = tmp_fmol -> first_atom;
2323 while (tmp_fat)
2324 {
2325 add_vdw = TRUE;
2326 for (j=0; j<i; j++)
2327 {
2328 if (g_strcmp0 (to_be_vdw[j], tmp_fat -> name) == 0)
2329 {
2330 add_vdw = FALSE;
2331 vdw_mlist[j] ++;
2332 vdw_mids[j][vdw_mlist[j]] = tmp_fmol -> id;
2333 vdw_aids[j][vdw_mlist[j]] = tmp_fat -> id;
2334 break;
2335 }
2336 }
2337 if (add_vdw)
2338 {
2339 to_be_vdw[j] = g_strdup_printf ("%s", tmp_fat -> name);
2340 vdw_mids[j][0] = tmp_fmol -> id;
2341 vdw_aids[j][0] = tmp_fat -> id;
2342 i ++;
2343 }
2344 tmp_fat = tmp_fat -> next;
2345 }
2346 tmp_fmol = tmp_fmol -> next;
2347 }
2348
2349 j = 0 ;
2350 tmp_fbody = NULL;
2351 for (k=0; k<i; k++)
2352 {
2353 vdw_na[0] = vdw_na[1] = vdw_mlist[k]+1;
2354 vdw_a[0] = vdw_a[1] = duplicate_int (vdw_na[0], vdw_aids[k]);
2355 vdw_ma[0] = vdw_ma[1] = duplicate_int (vdw_na[0], vdw_mids[k]);
2356 if (tmp_fbody == NULL)
2357 {
2358 if (init)
2359 {
2360 tmp_field -> first_body[0] = init_field_nth_body (j, 0, vdw_na, vdw_ma, vdw_a);
2361 tmp_fbody = tmp_field -> first_body[0];
2362 }
2363 else
2364 {
2365 comp_fbody = init_field_nth_body (j, 0, vdw_na, vdw_ma, vdw_a);
2366 tmp_fbody = comp_fbody;
2367 }
2368 }
2369 else
2370 {
2371 tmp_fbody -> next = init_field_nth_body (j, 0, vdw_na, vdw_ma, vdw_a);
2372 tmp_fbody -> next -> prev = tmp_fbody;
2373 tmp_fbody = tmp_fbody -> next;
2374 }
2375 j ++;
2376 }
2377 for (k=0; k<i-1; k++)
2378 {
2379 for (l=k+1; l<i; l++)
2380 {
2381 vdw_na[0] = vdw_mlist[k]+1;
2382 vdw_na[1] = vdw_mlist[l]+1;
2383 vdw_a[0] = duplicate_int (vdw_na[0], vdw_aids[k]);
2384 vdw_a[1] = duplicate_int (vdw_na[1], vdw_aids[l]);
2385 vdw_ma[0] = duplicate_int (vdw_na[0], vdw_mids[k]);
2386 vdw_ma[1] = duplicate_int (vdw_na[1], vdw_mids[l]);
2387 tmp_fbody -> next = init_field_nth_body (j, 0, vdw_na, vdw_ma, vdw_a);
2388 tmp_fbody -> next -> prev = tmp_fbody;
2389 tmp_fbody = tmp_fbody -> next;
2390 j ++;
2391 }
2392 }
2393 g_free (to_be_vdw);
2394 g_free (vdw_aids);
2395 g_free (vdw_mids);
2396 g_free (vdw_mlist);
2397 return j;
2398}
2399
2407void init_dlpoly_field (classical_field * new_field)
2408{
2409 // Preparing DLPOLY Control file
2410 new_field -> sys_opts = allocdouble (17);
2411 new_field -> io_opts = allocdouble (23);
2412 new_field -> ana_opts = allocdouble (17);
2413 new_field -> elec_opts = allocdouble (11);
2414 new_field -> vdw_opts = allocdouble (6);
2415 new_field -> met_opts = allocdouble (2);
2416 new_field -> equi_opts = allocdouble (17);
2417 new_field -> thermo_opts= allocdouble (10);
2418 new_field -> md_opts = allocdouble (20);
2419 new_field -> out_opts = allocdouble (31);
2420/*
2421 new_field -> io_opts[]
2422 */
2423 // Set job time control (yes/no)
2424 new_field -> io_opts[0] = 1.0;
2425 // Job time control (s):
2426 new_field -> io_opts[1] = 10000.0;
2427 // Set job closure time (yes/no)
2428 new_field -> io_opts[2] = 1.0;
2429 // Job closure time (s):
2430 new_field -> io_opts[3] = 100.0;
2431 // 4 (1/0) General I/O read interace
2432 // 5 General I/O method (0= MPI, 1= FORTRAN, 2= Traditional, 3= netCDF)
2433 // Reader count:
2434 new_field -> io_opts[6] = 12;
2435 // Batch size:
2436 new_field -> io_opts[7] = 2000000;
2437 // Buffer size:
2438 new_field -> io_opts[8] = 20000;
2439 // 9 (1/0) Error check on general read I/O
2440 // 10 (1/0) General I/O write interace
2441 // 11 General I/O method (0= MPI, 1= FORTRAN, 2= Traditional, 3= netCDF)
2442 // Real precision (0 = 32 bits, 1= 64 bits)
2443 new_field -> io_opts[12] = 1;
2444 // 13 Type (0= Sorted, 1= Unsorted)
2445 // Reader count:
2446 new_field -> io_opts[14] = 12;
2447 // Batch size:
2448 new_field -> io_opts[15] = 2000000;
2449 // Buffer size:
2450 new_field -> io_opts[16] = 20000;
2451 // 17 (1/0) Error check on general write I/O
2452 // 18 (1/0) Define random number generator seeds
2453 // n1 seed:
2454 new_field -> io_opts[19] = 1;
2455 // n2 seed:
2456 new_field -> io_opts[20] = 20;
2457 // n3 seed:
2458 new_field -> io_opts[21] = 300;
2459 // 22 (1/0) limit the number of processors in the z direction
2460
2461 /*
2462 new_field -> sys_opts[]
2463 */
2464 // Relative dieletric constant:
2465 new_field -> sys_opts[0] = 1.0;
2466 // 1 = allowed local variation of the system density (%)
2467 // 2 (1/0) Ingore the particle indices in CONFIG file
2468 // 3 (1/0) Ignore strict checks
2469 // 4 (1/0) Skip detailed topology reporting
2470 // 5 (1/0) Ignore center of mass removal
2471 // 6 (1/0) Set tolerance of the relaxed shell model
2472 // Tolerance for the relaxed shell model:
2473 new_field -> sys_opts[7] = 1.0;
2474 // 8 (1/0) Set the subscelling threshold density of particles per link cell
2475 // Subscelling threshold density of particles per link cell:
2476 new_field -> sys_opts[9] = 50.0;
2477 // 10 (1/0) Create an expanded version of the current model ...
2478 // ... on x:
2479 new_field -> sys_opts[11] = 1.0;
2480 // ... on y:
2481 new_field -> sys_opts[12] = 1.0;
2482 // .. on z:
2483 new_field -> sys_opts[13] = 1.0;
2484 // 14 (1/0) Restart a calculation
2485 // 15 (0= Continue current, 1= Start new without T reset, 2= Start new with T reset)
2486
2487 /*
2488 new_field -> equi_opts[]
2489 */
2490 new_field -> equi_opts[0] = 1.0;
2491 // Equilibrate for the first N steps, N=
2492 new_field -> equi_opts[1] = 20.0;
2493 // Scale temperature every N steps, N=
2494 new_field -> equi_opts[3] = 1.0;
2495 // Cap forces fmax=
2496 new_field -> equi_opts[5] = 1000.0;
2497 // Resample the momenta distribution during N steps N=
2498 new_field -> equi_opts[7] = 1.0;
2499 // 8 (1/0) = Minimize system configuration
2500 // 9 Minimize target (0=Force, 1=Energy, 2= Distance)
2501 // Minimization threshold:
2502 new_field -> equi_opts[10] = 50.0;
2503 // Minimization to be done every N steps, N=
2504 new_field -> equi_opts[11] = 1.0;
2505 // 12 (1/0) = Optimize system configuration from start
2506 // 13 Minimize target (0=Force, 1=Energy, 2= Distance)
2507 // Minimization threshold:
2508 new_field -> equi_opts[14] = 50.0;
2509 // 15 (1/0) = Perform zero temperature (10K) optimization
2510 // 16 (1/0) = Include equilibration data in statistics
2511
2512 /*
2513 new_field -> vdw_opts[]
2514 */
2515 // Vdw cutoff
2516 // 0 (1/0) Use VdW interactions
2517 new_field -> vdw_opts[0] = 1.0;
2518 // Vdw Cutoff radius
2519 new_field -> vdw_opts[1] = 12.0;
2520 // 2 (1/0) Enforce direct VdW calculation
2521 // 3 (1/0) Apply force shifting
2522 // 4 (1/0) Apply mixing rule
2523 // 5 Mixing rule (0= LB, 1= FH, 2= H, 3= HHG, 4= TT, 5= F)
2524
2525 /*
2526 new_field -> elec_opts[]
2527 */
2528 // 0 (1/0) Evaluate electrostatics
2529 new_field -> elec_opts[0] = 1.0;
2530 // Long range cutoff:
2531 new_field -> elec_opts[1] = 12.0;
2532 // 2 (1/0) Add padding to rcut
2533 // 3 Padding in angstrom
2534 // 4 (1/0) Use extended coulombic exclusion
2535 // 5 Electrostatic calculation method (0= Direct, ...)
2536 new_field -> elec_opts[5] = 2.0;
2537 // Ewald precision:
2538 new_field -> elec_opts[6] = 1e-20;
2539 // 7 to 9 method options
2540 // 10 Evaluation every N steps, N=
2541 tmp_field -> elec_opts[10] = 1.0;
2542
2543 /*
2544 new_field -> met_opts[]
2545 */
2546 // 0 (1/0) Enable direct calculation of metal interactions
2547 // 1 (1/0) Switch the default embeding fuctions
2548
2549 // Ensemble (0= NVE, 1= NVT, 2= NPT, 4= NST)
2550 // new_field -> ensemble
2551 // Thermostat = f(ensemble):
2552 // new_field -> thermostat
2553 /*
2554 new_field -> thermo_opts[]
2555 */
2556 // 0-5 ensemble + thermostat related options
2557 // 6 (1/0) Use pseudo thermal bath
2558 // 7 Thermostat type (0= L+D, 1= L, 2= G, 3= D)
2559 // Thickness:
2560 new_field -> thermo_opts[8] = 2.0;
2561 // 9 = Target temperature
2562
2563 /*
2564 new_field -> md_opts[]
2565 */
2566 // MD target temperature (K)
2567 new_field -> md_opts[0] = 300.0;
2568 // 1 Verlet integrator (0 = velocity, 1 = leapfrog)
2569 // Number of MD steps:
2570 new_field -> md_opts[2] = 100.0;
2571 // 3 Time step (0 = fixed, 1 = variable)
2572 // Time step delta(t)= (ps)
2573 new_field -> md_opts[4] = 0.001;
2574 // Target pressure (katms):
2575 new_field -> md_opts[5] = 0.001;
2576 // Variable time step, maximum time step:
2577 new_field -> md_opts[6] = 0.001;
2578 // Variable time step, maximum distance allowed:
2579 new_field -> md_opts[7] = 0.5;
2580 // Variable time step, minimun distance allowed:
2581 new_field -> md_opts[8] = 0.01;
2582 // Shake iterations limit, in cycles:
2583 new_field -> md_opts[9] = 250.0;
2584 // Shake tolerance:
2585 new_field -> md_opts[10] = 1e-6;
2586 // FIQA iterations limit, in cycles:
2587 new_field -> md_opts[11] = 100;
2588 // FIQA quaternions tolerance:
2589 new_field -> md_opts[12] = 1e-8;
2590 // 13 (1/0) Initiate impact on particles
2591 // Index if the particle to impact:
2592 new_field -> md_opts[14] = 1.0;
2593 // Time step of the impact:
2594 new_field -> md_opts[15] = 1.0;
2595 // Energy of the imapct (keV):
2596 new_field -> md_opts[16] = 1.0;
2597 // Direction of the impact (from center of mass) ...
2598 // ... on x:
2599 new_field -> md_opts[17] = 1.0;
2600 // ... on y:
2601 new_field -> md_opts[18] = 1.0;
2602 // ... on z:
2603 new_field -> md_opts[19] = 1.0;
2604
2605 /*
2606 new_field -> out_opts[]
2607 */
2608 // 0 (1/0) Print Defect file ...
2609 // 1 ... starting from step
2610 // ... every N steps, N=
2611 new_field -> out_opts[2] = 1.0;
2612 // ... with a site inteerstitial cutoff:
2613 new_field -> out_opts[3] = 0.75;
2614 // 4 (1/0) Print Displacement file ..
2615 // 5 ... starting from step
2616 // ... every N steps, N=
2617 new_field -> out_opts[6] = 1.0;
2618 // ... with a qualifying cutoff:
2619 new_field -> out_opts[7] = 0.15;
2620 // 8 (1/0) Print HISTORY file ...
2621 // 9 ... starting from step
2622 // ... every N steps, N=
2623 new_field -> out_opts[10] = 1.0;
2624 // 11 Data level (0, 1, 2)
2625 // 12 (1/0) Print Atom's MSD ...
2626 // 13 ... starting from step
2627 // ... every N steps, N=
2628 new_field -> out_opts[14] = 1.0;
2629 // 15 (1/0) Print system data ...
2630 // ... every N steps, N=
2631 new_field -> out_opts[16] = 1.0;
2632 // 17 (1/0) Accumulate statistics ...
2633 // ... every N steps, N=
2634 new_field -> out_opts[18] = 1.0;
2635 // 19 (1/0) Rolling average stack ...
2636 // ... to N steps, N=
2637 new_field -> out_opts[20] = 1.0;
2638 // 21 (1/0) RDF ...
2639 // ... every N steps, N=
2640 new_field -> out_opts[22] = 1.0;
2641 // ... with a bin size=
2642 new_field -> out_opts[23] = 0.05;
2643 // 24 (1/0) Overall VAF analysis ...
2644 // ... every N steps, N=
2645 new_field -> out_opts[25] = 1.0;
2646 // ... with a bin size=
2647 new_field -> out_opts[26] = 50.0;
2648 // 27 (1/0) Z density profile ...
2649 // ... every N steps, N=
2650 new_field -> out_opts[28] = 1.0;
2651 // 29 (1/0) ignoring time-av. VAFs
2652 // Dump restart information every N steps, N=
2653 new_field -> out_opts[30] = 1000.0;
2654
2655 /*
2656 new_field -> ana_opts[]
2657 */
2658 // 0 (1/0) Activate all intra-molecular PDF analysis ...
2659 // ... every N steps, N=
2660 new_field -> ana_opts[1] = 1.0;
2661 // ... number of delta r in:
2662 new_field -> ana_opts[2] = 100.0;
2663 // ... Cutoff:
2664 new_field -> ana_opts[3] = 2.0;
2665 // 4 (1/0) Activate bonds PDF analysis ...
2666 // ... every N steps, N=
2667 new_field -> ana_opts[5] = 1.0;
2668 // ... number of delta r in:
2669 new_field -> ana_opts[6] = 100.0;
2670 // ... Cutoff:
2671 new_field -> ana_opts[7] = 2.0;
2672 // 8 (1/0) Activate angles analysis ...
2673 // ... every N steps, N=
2674 new_field -> ana_opts[9] = 1.0;
2675 // ... number of delta angle in:
2676 new_field -> ana_opts[10] = 90.0;
2677 // 11 (1/0) Activate dihedrals analysis ...
2678 // ... every N steps, N=
2679 new_field -> ana_opts[12] = 1.0;
2680 // ... number of delta angle in:
2681 new_field -> ana_opts[13] = 90.0;
2682 // 14 (1/0) Activate inversions analysis ...
2683 // ... every N steps, N=
2684 new_field -> ana_opts[15] = 1.0;
2685 // ... number of delta angle in:
2686 new_field -> ana_opts[16] = 90.0;
2687
2688}
2689
2697void init_lammps_field (classical_field * new_field)
2698{
2699 // Preparing LAMMPS Control file
2700 new_field -> sys_opts = allocdouble (17);
2701 new_field -> io_opts = allocdouble (23);
2702 new_field -> ana_opts = allocdouble (17);
2703 new_field -> elec_opts = allocdouble (11);
2704 new_field -> vdw_opts = allocdouble (6);
2705 new_field -> met_opts = allocdouble (2);
2706 new_field -> equi_opts = allocdouble (17);
2707 new_field -> thermo_opts= allocdouble (10);
2708 new_field -> md_opts = allocdouble (20);
2709 new_field -> out_opts = allocdouble (31);
2710/*
2711 new_field -> io_opts[]
2712 */
2713 // Set job time control (yes/no)
2714 new_field -> io_opts[0] = 1.0;
2715 // Job time control (s):
2716 new_field -> io_opts[1] = 10000.0;
2717 // Set job closure time (yes/no)
2718 new_field -> io_opts[2] = 1.0;
2719 // Job closure time (s):
2720 new_field -> io_opts[3] = 100.0;
2721 // 4 (1/0) General I/O read interace
2722 // 5 General I/O method (0= MPI, 1= FORTRAN, 2= Traditional, 3= netCDF)
2723 // Reader count:
2724 new_field -> io_opts[6] = 12;
2725 // Batch size:
2726 new_field -> io_opts[7] = 2000000;
2727 // Buffer size:
2728 new_field -> io_opts[8] = 20000;
2729 // 9 (1/0) Error check on general read I/O
2730 // 10 (1/0) General I/O write interace
2731 // 11 General I/O method (0= MPI, 1= FORTRAN, 2= Traditional, 3= netCDF)
2732 // Real precision (0 = 32 bits, 1= 64 bits)
2733 new_field -> io_opts[12] = 1;
2734 // 13 Type (0= Sorted, 1= Unsorted)
2735 // Reader count:
2736 new_field -> io_opts[14] = 12;
2737 // Batch size:
2738 new_field -> io_opts[15] = 2000000;
2739 // Buffer size:
2740 new_field -> io_opts[16] = 20000;
2741 // 17 (1/0) Error check on general write I/O
2742 // 18 (1/0) Define random number generator seeds
2743 // n1 seed:
2744 new_field -> io_opts[19] = 1;
2745 // n2 seed:
2746 new_field -> io_opts[20] = 20;
2747 // n3 seed:
2748 new_field -> io_opts[21] = 300;
2749 // 22 (1/0) limit the number of processors in the z direction
2750
2751 /*
2752 new_field -> sys_opts[]
2753 */
2754 // Relative dieletric constant:
2755 new_field -> sys_opts[0] = 1.0;
2756 // 1 = allowed local variation of the system density (%)
2757 // 2 (1/0) Ingore the particle indices in CONFIG file
2758 // 3 (1/0) Ignore strict checks
2759 // 4 (1/0) Skip detailed topology reporting
2760 // 5 (1/0) Ignore center of mass removal
2761 // 6 (1/0) Set tolerance of the relaxed shell model
2762 // Tolerance for the relaxed shell model:
2763 new_field -> sys_opts[7] = 1.0;
2764 // 8 (1/0) Set the subscelling threshold density of particles per link cell
2765 // Subscelling threshold density of particles per link cell:
2766 new_field -> sys_opts[9] = 50.0;
2767 // 10 (1/0) Create an expanded version of the current model ...
2768 // ... on x:
2769 new_field -> sys_opts[11] = 1.0;
2770 // ... on y:
2771 new_field -> sys_opts[12] = 1.0;
2772 // .. on z:
2773 new_field -> sys_opts[13] = 1.0;
2774 // 14 (1/0) Restart a calculation
2775 // 15 (0= Continue current, 1= Start new without T reset, 2= Start new with T reset)
2776
2777 /*
2778 new_field -> equi_opts[]
2779 */
2780 new_field -> equi_opts[0] = 1.0;
2781 // Equilibrate for the first N steps, N=
2782 new_field -> equi_opts[1] = 20.0;
2783 // Scale temperature every N steps, N=
2784 new_field -> equi_opts[3] = 1.0;
2785 // Cap forces fmax=
2786 new_field -> equi_opts[5] = 1000.0;
2787 // Resample the momenta distribution during N steps N=
2788 new_field -> equi_opts[7] = 1.0;
2789 // 8 (1/0) = Minimize system configuration
2790 // 9 Minimize target (0=Force, 1=Energy, 2= Distance)
2791 // Minimization threshold:
2792 new_field -> equi_opts[10] = 50.0;
2793 // Minimization to be done every N steps, N=
2794 new_field -> equi_opts[11] = 1.0;
2795 // 12 (1/0) = Optimize system configuration from start
2796 // 13 Minimize target (0=Force, 1=Energy, 2= Distance)
2797 // Minimization threshold:
2798 new_field -> equi_opts[14] = 50.0;
2799 // 15 (1/0) = Perform zero temperature (10K) optimization
2800 // 16 (1/0) = Include equilibration data in statistics
2801
2802 /*
2803 new_field -> vdw_opts[]
2804 */
2805 // Vdw cutoff
2806 // 0 (1/0) Use VdW interactions
2807 new_field -> vdw_opts[0] = 1.0;
2808 // Vdw Cutoff radius
2809 new_field -> vdw_opts[1] = 12.0;
2810 // 2 (1/0) Enforce direct VdW calculation
2811 // 3 (1/0) Apply force shifting
2812 // 4 (1/0) Apply mixing rule
2813 // 5 Mixing rule (0= LB, 1= FH, 2= H, 3= HHG, 4= TT, 5= F)
2814
2815 /*
2816 new_field -> elec_opts[]
2817 */
2818 // 0 (1/0) Evaluate electrostatics
2819 new_field -> elec_opts[0] = 1.0;
2820 // Long range cutoff:
2821 new_field -> elec_opts[1] = 12.0;
2822 // 2 (1/0) Add padding to rcut
2823 // 3 Padding in angstrom
2824 // 4 (1/0) Use extended coulombic exclusion
2825 // 5 Electrostatic calculation method (0= Direct, ...)
2826 new_field -> elec_opts[5] = 2.0;
2827 // Ewald precision:
2828 new_field -> elec_opts[6] = 1e-20;
2829 // 7 to 9 method options
2830 // 10 Evaluation every N steps, N=
2831 tmp_field -> elec_opts[10] = 1.0;
2832
2833 /*
2834 new_field -> met_opts[]
2835 */
2836 // 0 (1/0) Enable direct calculation of metal interactions
2837 // 1 (1/0) Switch the default embeding fuctions
2838
2839 // Ensemble (0= NVE, 1= NVT, 2= NPT, 4= NST)
2840 // new_field -> ensemble
2841 // Thermostat = f(ensemble):
2842 // new_field -> thermostat
2843 /*
2844 new_field -> thermo_opts[]
2845 */
2846 // 0-5 ensemble + thermostat related options
2847 // 6 (1/0) Use pseudo thermal bath
2848 // 7 Thermostat type (0= L+D, 1= L, 2= G, 3= D)
2849 // Thickness:
2850 new_field -> thermo_opts[8] = 2.0;
2851 // 9 = Target temperature
2852
2853 /*
2854 new_field -> md_opts[]
2855 */
2856 // MD target temperature (K)
2857 new_field -> md_opts[0] = 300.0;
2858 // 1 Verlet integrator (0 = velocity, 1 = leapfrog)
2859 // Number of MD steps:
2860 new_field -> md_opts[2] = 100.0;
2861 // 3 Time step (0 = fixed, 1 = variable)
2862 // Time step delta(t)= (ps)
2863 new_field -> md_opts[4] = 0.001;
2864 // Target pressure (katms):
2865 new_field -> md_opts[5] = 0.001;
2866 // Variable time step, maximum time step:
2867 new_field -> md_opts[6] = 0.001;
2868 // Variable time step, maximum distance allowed:
2869 new_field -> md_opts[7] = 0.5;
2870 // Variable time step, minimun distance allowed:
2871 new_field -> md_opts[8] = 0.01;
2872 // Shake iterations limit, in cycles:
2873 new_field -> md_opts[9] = 250.0;
2874 // Shake tolerance:
2875 new_field -> md_opts[10] = 1e-6;
2876 // FIQA iterations limit, in cycles:
2877 new_field -> md_opts[11] = 100;
2878 // FIQA quaternions tolerance:
2879 new_field -> md_opts[12] = 1e-8;
2880 // 13 (1/0) Initiate impact on particles
2881 // Index if the particle to impact:
2882 new_field -> md_opts[14] = 1.0;
2883 // Time step of the impact:
2884 new_field -> md_opts[15] = 1.0;
2885 // Energy of the imapct (keV):
2886 new_field -> md_opts[16] = 1.0;
2887 // Direction of the impact (from center of mass) ...
2888 // ... on x:
2889 new_field -> md_opts[17] = 1.0;
2890 // ... on y:
2891 new_field -> md_opts[18] = 1.0;
2892 // ... on z:
2893 new_field -> md_opts[19] = 1.0;
2894
2895 /*
2896 new_field -> out_opts[]
2897 */
2898 // 0 (1/0) Print Defect file ...
2899 // 1 ... starting from step
2900 // ... every N steps, N=
2901 new_field -> out_opts[2] = 1.0;
2902 // ... with a site inteerstitial cutoff:
2903 new_field -> out_opts[3] = 0.75;
2904 // 4 (1/0) Print Displacement file ..
2905 // 5 ... starting from step
2906 // ... every N steps, N=
2907 new_field -> out_opts[6] = 1.0;
2908 // ... with a qualifying cutoff:
2909 new_field -> out_opts[7] = 0.15;
2910 // 8 (1/0) Print HISTORY file ...
2911 // 9 ... starting from step
2912 // ... every N steps, N=
2913 new_field -> out_opts[10] = 1.0;
2914 // 11 Data level (0, 1, 2)
2915 // 12 (1/0) Print Atom's MSD ...
2916 // 13 ... starting from step
2917 // ... every N steps, N=
2918 new_field -> out_opts[14] = 1.0;
2919 // 15 (1/0) Print system data ...
2920 // ... every N steps, N=
2921 new_field -> out_opts[16] = 1.0;
2922 // 17 (1/0) Accumulate statistics ...
2923 // ... every N steps, N=
2924 new_field -> out_opts[18] = 1.0;
2925 // 19 (1/0) Rolling average stack ...
2926 // ... to N steps, N=
2927 new_field -> out_opts[20] = 1.0;
2928 // 21 (1/0) RDF ...
2929 // ... every N steps, N=
2930 new_field -> out_opts[22] = 1.0;
2931 // ... with a bin size=
2932 new_field -> out_opts[23] = 0.05;
2933 // 24 (1/0) Overall VAF analysis ...
2934 // ... every N steps, N=
2935 new_field -> out_opts[25] = 1.0;
2936 // ... with a bin size=
2937 new_field -> out_opts[26] = 50.0;
2938 // 27 (1/0) Z density profile ...
2939 // ... every N steps, N=
2940 new_field -> out_opts[28] = 1.0;
2941 // 29 (1/0) ignoring time-av. VAFs
2942 // Dump restart information every N steps, N=
2943 new_field -> out_opts[30] = 1000.0;
2944
2945 /*
2946 new_field -> ana_opts[]
2947 */
2948 // 0 (1/0) Activate all intra-molecular PDF analysis ...
2949 // ... every N steps, N=
2950 new_field -> ana_opts[1] = 1.0;
2951 // ... number of delta r in:
2952 new_field -> ana_opts[2] = 100.0;
2953 // ... Cutoff:
2954 new_field -> ana_opts[3] = 2.0;
2955 // 4 (1/0) Activate bonds PDF analysis ...
2956 // ... every N steps, N=
2957 new_field -> ana_opts[5] = 1.0;
2958 // ... number of delta r in:
2959 new_field -> ana_opts[6] = 100.0;
2960 // ... Cutoff:
2961 new_field -> ana_opts[7] = 2.0;
2962 // 8 (1/0) Activate angles analysis ...
2963 // ... every N steps, N=
2964 new_field -> ana_opts[9] = 1.0;
2965 // ... number of delta angle in:
2966 new_field -> ana_opts[10] = 90.0;
2967 // 11 (1/0) Activate dihedrals analysis ...
2968 // ... every N steps, N=
2969 new_field -> ana_opts[12] = 1.0;
2970 // ... number of delta angle in:
2971 new_field -> ana_opts[13] = 90.0;
2972 // 14 (1/0) Activate inversions analysis ...
2973 // ... every N steps, N=
2974 new_field -> ana_opts[15] = 1.0;
2975 // ... number of delta angle in:
2976 new_field -> ana_opts[16] = 90.0;
2977}
2978
2986classical_field * create_force_field_data_structure (int ai)
2987{
2988 int i, j;
2989
2990 classical_field * new_field;
2991 new_field = g_malloc0 (sizeof*new_field);
2992 tmp_field = new_field;
2993 new_field -> energy_unit = 1;
2994 new_field -> atom_init = ai;
2995 switch (activef)
2996 {
2997 case 0:
2998 init_dlpoly_field (new_field);
2999 break;
3000 case 1:
3001 init_lammps_field (new_field);
3002 break;
3003 }
3004 new_field -> molecules = tmp_proj -> modelfc -> mol_by_step[0];
3005 new_field -> first_molecule = g_malloc0 (sizeof*new_field -> first_molecule);
3006 tmp_fmol = new_field -> first_molecule;
3007 j = 0;
3008#ifdef DEBUG
3009 g_debug ("Total num of molecules: %d", tmp_proj -> modelfc -> mol_by_step[0]);
3010#endif // DEBUG
3011 for (i=0; i < tmp_proj -> modelfc -> mol_by_step[0]; i++)
3012 {
3013 tmp_mol = & tmp_proj -> modelfc -> mols[0][i];
3014#ifdef DEBUG
3015 g_debug ("Init FORCE FIELD for: ");
3016 g_debug (" tmp_mol -> id= %d", tmp_mol -> id);
3017 g_debug (" tmp_mol -> multi= %d", tmp_mol -> multiplicity);
3018 g_debug (" tmp_mol -> natoms= %d", tmp_mol -> natoms);
3019#endif // DEBUG
3020 init_field_molecule (i);
3021 j += tmp_fmol -> atoms;
3022 if (i < tmp_proj -> modelfc -> mol_by_step[0]-1)
3023 {
3024 tmp_fmol -> next = g_malloc0 (sizeof*tmp_fmol -> next);
3025 tmp_fmol = tmp_fmol -> next;
3026 }
3027 }
3028
3029 for (i=0; i<5; i++)
3030 {
3031 tmp_field -> nbody[i] = 0;
3032 tmp_field -> first_body[i] = NULL;
3033 }
3034 // Tersoff cross terms
3035 tmp_field -> cross = NULL;
3036
3037
3038 // vdW start
3039 tmp_field -> nbody[0] = init_vdw (TRUE);
3040 // vdW end
3041
3042 // External field(s) start
3043 new_field -> extern_fields = 0;
3044 new_field -> first_external = NULL;
3045 // External field(s) end
3046
3047 // Non-bonded end
3048 return new_field;
3049}
limit
float limit[2]
Definition atom_edit.c:59
mol
insertion_menu mol[]
Definition w_library.c:193
search_type
int search_type
Definition calc_menu.c:91
chemistry
integer(kind=c_int) function chemistry()
Definition chemistry.F90:22
label
void label(cairo_t *cr, double val, int axe, int p, project *this_proj)
draw axis label
Definition labels.c:56
atoms
int atoms[NUM_STYLES][2]
Definition d_selection.c:361
get_active_atom
field_atom * get_active_atom(int a, int b)
retrieve field atom
Definition dlp_active.c:145
comp_fbody
field_nth_body * comp_fbody
Definition dlp_comp.c:44
sys_opts
gchar * sys_opts[10]
Definition dlp_control.c:2410
duplicate_field_struct_list
field_struct * duplicate_field_struct_list(field_struct *list_str, gboolean init)
create copy of list of field structural element(s)
Definition dlp_copy.c:325
atoms_id
int * atoms_id
Definition force_fields.c:238
cross
double *** cross
Definition dlp_edit.c:418
fvalues
int fvalues[2][15][21]
Definition dlp_field.c:255
tmp_fbody
field_nth_body * tmp_fbody
Definition dlp_field.c:965
tmp_proj
project * tmp_proj
Definition dlp_field.c:953
body_at
int body_at(int b)
find the number of atom(s) in a non bonded interaction
Definition dlp_field.c:1022
tmp_fat
field_atom * tmp_fat
Definition dlp_field.c:957
struct_id
int struct_id(int f)
number of atoms in a structural element
Definition dlp_field.c:999
tmp_fstr
field_struct * tmp_fstr
Definition dlp_field.c:964
tmp_fmol
field_molecule * tmp_fmol
Definition dlp_field.c:955
atom_init
int atom_init
Definition dlp_field.c:950
tmp_fdt
field_atom * tmp_fdt
Definition dlp_field.c:957
tmp_fct
field_atom * tmp_fct
Definition dlp_field.c:957
tmp_fbt
field_atom * tmp_fbt
Definition dlp_field.c:957
tmp_coord
coord_info * tmp_coord
Definition dlp_field.c:954
tmp_mol
molecule * tmp_mol
Definition dlp_field.c:956
tmp_field
classical_field * tmp_field
Definition dlp_field.c:951
dlp_field.h
Variable declarations for the creation of the DL_POLY input file(s)
dihedral_3d
angle dihedral_3d(cell_info *cell, int mdstep, atom *at, atom *bt, atom *ct, atom *dt)
dihedral between atom a, b, c and d in 3D
Definition ogl_utils.c:204
inversion_3d
angle inversion_3d(cell_info *cell, int mdstep, atom *at, atom *bt, atom *ct, atom *dt)
inversion angle between atom a, b, c and d in 3D
Definition ogl_utils.c:240
angle_3d
angle angle_3d(cell_info *cell, int mdstep, atom *at, atom *bt, atom *ct)
angle between atom a, b and c in 3D
Definition ogl_utils.c:179
TOTAL_COORD_AND_SPEC
@ TOTAL_COORD_AND_SPEC
Definition dlp_field.h:146
SPEC_ONLY
@ SPEC_ONLY
Definition dlp_field.h:144
PARTIAL_COORD_AND_SPEC
@ PARTIAL_COORD_AND_SPEC
Definition dlp_field.h:148
OTHER
@ OTHER
Definition dlp_field.h:149
distance_3d
distance distance_3d(cell_info *cell, int mdstep, atom *at, atom *bt)
distance between atom a and b in 3D
Definition ogl_utils.c:81
angles_from_bonds
int angles_from_bonds(int n, field_atom *at, field_atom *bt, field_atom *ct)
find, and initialize, angles using bonds
Definition dlp_init.c:1129
init_field_shell
field_shell * init_field_shell(int id, int ia, int ib)
initialize field core shell interaction
Definition dlp_init.c:394
init_field_tethered
field_tethered * init_field_tethered(int id, int num)
intialize new field tethered potential
Definition dlp_init.c:510
id_atom
int id_atom
Definition dlp_init.c:120
clean_field_struct_list
void clean_field_struct_list(field_struct *stru)
clean force field structural property(ies)
Definition dlp_init.c:2146
find_neighbor_loop
int find_neighbor_loop(int frag, int aid, int fai, field_atom *fat, field_neighbor *ngmb, int sid, gboolean save_it)
find / set up field atom data for a fragment
Definition dlp_init.c:1591
init_ngb
field_neighbor * init_ngb
Definition dlp_init.c:1474
init_dlpoly_field
void init_dlpoly_field(classical_field *new_field)
initialize DL-POLY classical force field
Definition dlp_init.c:2407
init_field_struct
field_struct * init_field_struct(int st, int ai, int an, int *aid)
initialize field molecule new structural element
Definition dlp_init.c:704
duplicate_int
int * duplicate_int(int num, int *old_val)
copy a list of int
Definition global.c:572
prepare_field_atom
int prepare_field_atom(int i, int j, int k, int l, int m)
initialiaze new type of field atom
Definition dlp_init.c:814
test_for_angles
int test_for_angles(field_atom *at, field_atom *bt, field_atom *ct)
search for angle(s) between these field atoms
Definition dlp_init.c:1081
assigned
int assigned
Definition dlp_init.c:1471
val
float val
Definition dlp_init.c:117
get_struct_id_from_atom_id
int get_struct_id_from_atom_id(int ids, int *aid)
retrieve field structural element id from a list of atom id
Definition dlp_init.c:559
get_atom_id_in_fragment_from_model_id
int get_atom_id_in_fragment_from_model_id(int frag, int at)
retrieve field atom id in fragment from model id
Definition dlp_init.c:197
id_n_fold_atoms_in_fragment
gboolean id_n_fold_atoms_in_fragment(int frag, int limit, int num_ngb, int search_type, int init, field_neighbor *ngma_init)
find atom id in fragment using the neighbor list
Definition dlp_init.c:1638
get_init_neighbor
field_neighbor * get_init_neighbor(int a)
retrieve neighbor data structure for atom a
Definition dlp_init.c:1483
print_all_field_struct
void print_all_field_struct(field_molecule *mol, int str)
print all field structural element(s)
Definition dlp_print.c:195
set_atom_id
int set_atom_id(field_atom *at, int c, int p, int id)
set new fragment/molecule atom parameters
Definition dlp_init.c:224
coord_sphere_multiplicity
int coord_sphere_multiplicity(atom *at, int id, gboolean set_atom_id)
get coord sphere multiplicity in the model
Definition dlp_init.c:1420
get_weight
void get_weight(int seq)
compute field weigth for the atom(s) in the model
Definition dlp_init.c:1832
init_all_impropers_inversions
void init_all_impropers_inversions(int stru)
initialize all impropers and inversions
Definition dlp_init.c:1343
set_field_atom_name
gchar * set_field_atom_name(field_atom *ato, field_molecule *mol)
get name string for field atom
Definition dlp_init.c:247
init_field_prop
field_prop * init_field_prop(int ti, int key, gboolean show, gboolean use)
initialize new field molecule structural property
Definition dlp_init.c:535
get_field_atom_id_from_model_id
int get_field_atom_id_from_model_id(field_molecule *fmol, int at)
retrieve field atom id using atom id in model
Definition dlp_init.c:151
init_field_molecule
void init_field_molecule(int i)
initialiaze of field molecule data structurre
Definition dlp_init.c:2239
multi
int multi
Definition dlp_init.c:121
are_in_bond
gboolean are_in_bond(atom ato, int at)
is at in ato neighbors ?
Definition dlp_init.c:930
was_not_created_struct
gboolean was_not_created_struct(int ids, int num, int *aid)
was this structural element already created ?
Definition dlp_init.c:603
astr
int * astr
Definition dlp_init.c:116
a_multi
int a_multi
Definition dlp_init.c:122
test_for_bonds
int test_for_bonds(field_atom *at, field_atom *bt)
search for bond(s) between 2 field atoms
Definition dlp_init.c:948
get_fragment_atom_id_from_model_id
int get_fragment_atom_id_from_model_id(field_molecule *fmol, int at)
retrieve fragment atom id using atom id in model
Definition dlp_init.c:174
i_weight
int * i_weight
Definition dlp_init.c:1472
bonds_between_atoms
int bonds_between_atoms(int n, field_atom *at, field_atom *bt, int a, int b)
search for bond(s) between two field atoms
Definition dlp_init.c:1015
setup_atomic_weight
int setup_atomic_weight(int seq)
setup atom(s) field weight
Definition dlp_init.c:1877
init_all_dihedrals
void init_all_dihedrals()
find, and intialiaze, all dihedral(s) using bonds
Definition dlp_init.c:1259
dihedrals_from_angles
int dihedrals_from_angles(int n, field_atom *at, field_atom *bt, field_atom *ct, field_atom *dt)
find, and initialize, dihedrals using angles
Definition dlp_init.c:1240
impropers_inversion
int impropers_inversion(int n, int stru, int at, int bt, int ct, int dt, int a, int b, int c, int d)
prepare new improper / inversion
Definition dlp_init.c:1301
init_field_pmf
field_pmf * init_field_pmf(int id, int num[2], int *list[2], float *w[2])
initialize new field mean force potential
Definition dlp_init.c:445
duplicate_float
float * duplicate_float(int num, float *old_val)
copy a list of float
Definition global.c:604
test_for_dihedrals
int test_for_dihedrals(field_atom *at, field_atom *bt, field_atom *ct, field_atom *dt)
search for dihedral(s) between these field atoms
Definition dlp_init.c:1181
create_force_field_data_structure
classical_field * create_force_field_data_structure(int ai)
initialize classical force field
Definition dlp_init.c:2986
n_weight
int * n_weight
Definition dlp_init.c:1473
get_position_in_field_atom_from_model_id
int get_position_in_field_atom_from_model_id(int fat, int at)
retrieve atom position id in field atom id using atom id in model
Definition dlp_init.c:132
prepare_field_struct
int prepare_field_struct(int ids, int sid, int yes_no_num, int *aid)
prepare the creation of a field structural element
Definition dlp_init.c:981
init_field_external
field_external * init_field_external(int bi)
intialize new field external potential
Definition dlp_init.c:790
find_atom_id_in_field_molecule
void find_atom_id_in_field_molecule()
find field atom id in field molecule (fragment multiplicity > 1)
Definition dlp_init.c:1972
are_neighbors
gboolean are_neighbors(field_neighbor *ngb, int at)
test if 2 atoms are neighbors
Definition dlp_init.c:912
init_all_atoms
void init_all_atoms(int i)
initialize all field atom(s)
Definition dlp_init.c:838
id_atoms_in_fragment
gboolean id_atoms_in_fragment(int frag, int seed)
find the atom field id in the fragment
Definition dlp_init.c:1789
in_bond
gboolean in_bond(int at, int bd[2])
is atom at in bond bd
Definition dlp_init.c:892
init_field_constraint
field_constraint * init_field_constraint(int id, int ia, int ib)
initialize field bond constraint
Definition dlp_init.c:419
init_all_angles
void init_all_angles()
find, and intialiaze, all angle(s) using bonds
Definition dlp_init.c:1145
init_field_nth_body
field_nth_body * init_field_nth_body(int bi, int bd, int *na, int **ma, int **ba)
intialize new field non bonded interaction
Definition dlp_init.c:740
init_field_atom
field_atom * init_field_atom(int id, int type, int at, int nat, int coo, int *list)
intialize a new type of field atom
Definition dlp_init.c:304
init_all_bonds
void init_all_bonds()
find, and initialize all bond(s)
Definition dlp_init.c:1031
init_field_rigid
field_rigid * init_field_rigid(int id, int num, int *list)
initialize new field rigid constraint
Definition dlp_init.c:487
init_lammps_field
void init_lammps_field(classical_field *new_field)
initialize a LAMMPS classical force field
Definition dlp_init.c:2697
init_all_field_struct
void init_all_field_struct(gboolean init)
initialize the force field structural properties for a field molecule
Definition dlp_init.c:2170
is_this_numbering_possible_for_this_atom
gboolean is_this_numbering_possible_for_this_atom(int frag, field_neighbor *ngma, int atom)
check if the id in fragment for this atom is possible or not
Definition dlp_init.c:1550
is_numbering_possible
gboolean is_numbering_possible(int frag)
compare if the atom numbering of two fragment(s) are identical, and it should
Definition dlp_init.c:1502
init_vdw
int init_vdw(gboolean init)
initialize the VdW section of a classical force field
Definition dlp_init.c:2299
allocdint
int ** allocdint(int xal, int yal)
allocate an int ** pointer
Definition global.c:342
allocbool
gboolean * allocbool(int val)
allocate a gboolean * pointer
Definition global.c:266
activef
int activef
Definition global.c:161
allocdouble
double * allocdouble(int val)
allocate a double * pointer
Definition global.c:471
allocint
int * allocint(int val)
allocate an int * pointer
Definition global.c:326
allocfloat
float * allocfloat(int val)
allocate a float * pointer
Definition global.c:410
CHEM_M
#define CHEM_M
Definition global.h:270
exact_name
gchar * exact_name(gchar *name)
short cut to print string without spaces
Definition interface.c:370
interface.h
Messaging function declarations.
molecules
integer(kind=c_int) function molecules(frag_and_mol, allbonds)
Definition molecules.F90:155
angle::angle
double angle
Definition glwin.h:115
atom
Definition global.h:839
atom::numv
int numv
Definition global.h:846
atom::vois
int * vois
Definition global.h:847
classical_field
Definition global.h:666
distance::length
double length
Definition glwin.h:123
field_atom
Definition dlp_field.h:168
field_constraint
Definition dlp_field.h:210
field_external
Definition dlp_field.h:299
field_molecule
Definition dlp_field.h:321
field_neighbor
Definition dlp_field.h:311
field_nth_body
Definition dlp_field.h:281
field_pmf
Definition dlp_field.h:223
field_prop
Definition dlp_field.h:154
field_rigid
Definition dlp_field.h:238
field_shell
Definition dlp_field.h:192
field_struct
Definition dlp_field.h:264
field_tethered
Definition dlp_field.h:250
project::coord
coord_info * coord
Definition global.h:908
b
int b
Definition tab-1.c:95
c
int c
Definition tab-1.c:95
d
int d
Definition tab-1.c:95
a
int a
Definition tab-1.c:95
res
GtkWidget * res[2]
Definition w_encode.c:212