[ViewVC] Diff of: freemol/trunk/src/mengine/src/pcm7.c

#	Line 3 \| Line 3
3		#include "pcwin.h"
4		#include "pcmod.h"
5
6	+	#include "attached.h"
7		#include "energies.h"
7	–	#include "angles.h"
8	–	#include "rings.h"
9	–	#include "torsions.h"
8		#include "bonds_ff.h"
9		#include "derivs.h"
12	–	#include "hess.h"
10		#include "pot.h"
14	–	#include "gmmx.h"
11		#include "utility.h"
12		#include "solv.h"
13	<	#include "field.h"
18	<
19	<	EXTERN struct t_optimize {
20	<	int param_avail, converge;
21	<	float initial_energy, final_energy, initial_heat, final_heat;
22	<	} optimize_data;
23	<	EXTERN struct t_units {
24	<	double bndunit, cbnd, qbnd;
25	<	double angunit, cang, qang, pang, sang, aaunit;
26	<	double stbnunit, ureyunit, torsunit, storunit, v14scale;
27	<	double aterm, bterm, cterm, dielec, chgscale;
28	<	} units;
13	>	#include "dipmom.h"
14
30	–	struct t_cbonds {
31	–	int nbnd, icbonds[MAXCBND][3], lpd[MAXCBND], ia1[MAXCBND],ia2[MAXCBND];
32	–	double vrad[MAXCBND],veps[MAXCBND];
33	–	} cbonds;
34	–	EXTERN struct t_high_coord {
35	–	int ncoord, i13[400][3];
36	–	} high_coord;
37	–
15		EXTERN struct t_minim_control {
16		int type, method, field, added_const;
17		char added_path[256],added_name[256];
18		} minim_control;
42	–
43	–	struct t_dipolemom {
44	–	double total, xdipole, ydipole, zdipole;
45	–	} dipolemom;
19
20		EXTERN struct t_minim_values {
21		int iprint, ndc, nconst;
#	Line 59 \| Line 32
32
33		FILE * fopen_path ( char * , char * , char * ) ;
34		void message_alert(char ,char );
62	–	void get_added_const(void);
63	–	void assign_gast_charge();
64	–	void lattice(void);
65	–	void bounds(void);
35		void pcm7_min(void);
67	–	int iscoord_bond(int, int);
68	–	void get_coordbonds(void);
36		int isbond(int, int);
70	–	int isangle(int,int);
71	–	int ishbond(int, int);
72	–	int istorsion(int, int);
73	–	void get_hbond(void);
74	–	void mark_hbond(void);
75	–	void reset_atom_data(void);
37		int initialise(void);
38		void get_bonds(void);
78	–	void get_angles(void);
39		void get_torsions(void);
40		void get_rings(void);
41		void set_field(void);
42		double energy(void);
43		double print_energy(void);
44		void pcmfout(int);
85	–	void orbital(void);
86	–	void free_orbit(void);
45		void set_active(void);
46	<	void estrtor(void);
47	<	void echarge(void);
48	<	void ewald(void);
49	<	void eurey(void);
50	<	void ebuck(void);
51	<	void ebuckmm3(void);
46	>	int setup_calculation(void);
47	>	void end_calculation(void);
48	>
49	>	void egeom1(void);
50	>	void esolv1(void);
51	>	void egeom(void);
52	>	void esolv(void);
53	>
54	>	void ehal(void);
55	>	void ehal1(void);
56	>	void ebufcharge(void);
57	>	void ebufcharge1(void);
58	>	void eopbend_wilson(void);
59	>	void eopbend_wilson1(void);
60	>
61		void estrbnd(void);
95	–	void edipole(void);
62		void etorsion(void);
63		void eangang(void);
64		void eangle(void);
65		void ebond(void);
66	<	void eopbend(void);
101	<	void ehbond(void);
102	<	void eimprop(void);
103	<	void eimptors(void);
104	<	void elj(void);
105	<	void ecoord_bond(void);
106	<	void elj_qq(void);
107	<	void ebuck_charge(void);
108	<	void ehigh_coord(void);
109	<
110	<	void estrtor1(void);
111	<	void echarge1(void);
112	<	void ebuck1(void);
113	<	void ebuckmm31(void);
66	>
67		void estrbnd1(void);
115	–	void edipole1(void);
68		void etorsion1(void);
117	–	void eangang1(void);
69		void eangle1(void);
70		void ebond1(void);
120	–	void eopbend1(void);
121	–	void ehbond1(void);
122	–	void eimprop1(void);
123	–	void eimptors1(void);
124	–	void elj1(void);
125	–	void ecoord_bond1(void);
126	–	void ewald1(void);
127	–	void eurey1(void);
128	–	void elj_qq1(void);
129	–	void ebuck_charge1(void);
130	–	void ehigh_coord1(void);
71
72		void kangle(void);
133	–	void korbit(void);
73		void ktorsion(void);
135	–	void kdipole(void);
74		void kstrbnd(void);
75		void kcharge(void);
76		void ksolv(void);
139	–	void kangang(void);
77		void kopbend(void);
141	–	void kstrtor(void);
142	–	void khbond(void);
143	–	void piseq(int, int);
144	–	void kimprop(void);
145	–	void kimptors(void);
146	–	void kcoord_bonds(void);
78		void kvdw(void);
79		int kbond(void);
149	–	void kewald(void);
150	–	void kurey(void);
80
81		void get_memory(void);
82		void free_memory(void);
83		void gradient(void);
84		void minimize(void);
156	–	void dynamics(void);
85		void read_datafiles(char *);
86	+	void get_added_const(void);
87		void attach(void);
159	–	void eheat(void);
160	–	void pirite(void);
161	–	void charge_dipole(void);
162	–	void dipole_dipole(void);
163	–	void piden(void);
88		void hdel(int);
89		void set_atomtypes(int);
90		void type(void);
91		void generate_bonds(void);
92		void zero_data(void);
169	–	void ehal(void);
170	–	void ehal1(void);
171	–	void ebufcharge(void);
172	–	void ebufcharge1(void);
173	–	void eopbend_wilson(void);
174	–	void eopbend_wilson1(void);
175	–	void hbondreset(void);
176	–	void vibrate(void);
177	–	void pireset(void);
178	–	void adjust_mmfftypes(void);
179	–	void build_neighbor_list(int);
180	–	int setup_calculation(void);
181	–	void end_calculation(void);
182	–	void egeom1(void);
183	–	void esolv1(void);
184	–	void egeom(void);
185	–	void esolv(void);
93		void xlogp(float *);
94
188	–	void mmxsub(int ia)
189	–	{
190	–	int nret;
191	–
192	–	minim_control.type = ia; // 0 to minimize, 1 for single point
193	–	nret = setup_calculation();
194	–	if (nret == FALSE)
195	–	{
196	–	end_calculation();
197	–	return;
198	–	}
199	–
200	–	pcm7_min();
201	–	end_calculation();
202	–	}
203	–	// ================================================================
204	–	void pcm7_min()
205	–	{
206	–	int i, print;
207	–	double etot;
208	–
209	–	print = FALSE;
210	–	if (minim_values.iprint == TRUE)
211	–	{
212	–	print = TRUE;
213	–	minim_values.iprint = FALSE;
214	–	}
215	–
216	–	optimize_data.initial_energy = energies.total;
217	–
218	–	if (minim_control.type == 1) // single point calculation
219	–	return;
220	–
221	–	minimize();
222	–
223	–	if (print == TRUE)
224	–	minim_values.iprint = TRUE;
225	–
226	–	if (minim_values.iprint == TRUE)
227	–	etot = print_energy();
228	–	else
229	–	etot = energy();
230	–
231	–	optimize_data.final_energy = energies.total;
232	–
233	–	// compute dipole moment here
234	–	dipolemom.xdipole = 0.0;
235	–	dipolemom.ydipole = 0.0;
236	–	dipolemom.zdipole = 0.0;
237	–	if (pot.use_charge \|\| pot.use_bufcharge) charge_dipole();
238	–	if (pot.use_dipole) dipole_dipole();
239	–	dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
240	–	dipolemom.ydipole*dipolemom.ydipole +
241	–	dipolemom.zdipole*dipolemom.zdipole);
242	–	}
95		// ==================================================================
96		int setup_calculation()
97		{
98	<	int nRet;
98	>	int i,j,nRet;
99		char string[30];
100		double etot;
101		if (minim_control.field == MM3)
#	Line 277 \| Line 129
129		generate_bonds();
130		if (minim_control.added_const)
131		get_added_const();
280	–	} else if (minim_control.field == AMBER)
281	–	{
282	–	pot.use_hbond = FALSE;
283	–	pot.use_picalc = FALSE;
284	–	hbondreset();
285	–	pireset();
286	–	hdel(2);
287	–	set_field();
288	–	generate_bonds();
289	–	if (minim_control.added_const)
290	–	get_added_const();
291	–	} else if (minim_control.field == OPLSAA)
292	–	{
293	–	pot.use_hbond = FALSE;
294	–	pot.use_picalc = FALSE;
295	–	hbondreset();
296	–	pireset();
297	–	hdel(1);
298	–	set_field();
299	–	generate_bonds();
300	–	if (minim_control.added_const)
301	–	get_added_const();
132		} else
133		{
134		if (minim_control.added_const)
#	Line 321 \| Line 151
151		get_bonds();
152		get_angles();
153		get_torsions();
324	–	get_coordbonds();
325	–
326	–	if (high_coord.ncoord > 0)
327	–	pot.use_highcoord = TRUE;
154
155		attach();
330	–
331	–
156		get_rings();
157		// need allene
158
159		// set active atoms
160		set_active();
161	<
161	>	// setup non_bonded list of atoms to skip
162	>	for (i=1; i <= natom; i++)
163	>	{
164	>	for (j=0; j < MAXIAT; j++)
165	>	if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
166	>	{
167	>	skip[i][atom[i].iat[j]] = i;
168	>	skip[atom[i].iat[j]][i] = atom[i].iat[j];
169	>	}
170	>	for (j=0; j < attached.n13[i]; j++)
171	>	{
172	>	skip[i][attached.i13[j][i]] = i;
173	>	skip[attached.i13[j][i]][i] = attached.i13[j][i];
174	>	}
175	>	for(j=0; j < attached.n14[i]; j++)
176	>	{
177	>	skip[i][attached.i14[j][i]] = -i;
178	>	skip[attached.i14[j][i]][i] = -attached.i14[j][i];
179	>	}
180	>	}
181		/* assign local geometry potential function parameters */
182		Missing_constants = FALSE;
183		// errfile = fopen_path(pcwindir,"pcmod.err","w");
#	Line 342 \| Line 185
185		if (nRet == FALSE) // use end_calc to free memory
186		{
187		// message_alert("Parameters missing. See pcmod.err for information","Error");
345	–	optimize_data.param_avail = 1;
188		// fclose(errfile);
189		energies.total = 9999.;
190		return FALSE;
#	Line 359 \| Line 201
201		{
202		// message_alert("Parameters missing. See pcmod.err for information","Error");
203		// fclose(errfile);
362	–	optimize_data.param_avail = 1;
204		energies.total = 9999.0;
205		return (FALSE);
206		} else
#	Line 401 \| Line 242
242		energies.eurey = 0.0;
243		energies.esolv = 0.0;
244		energies.egeom =0.0;
404	–
405	–	virial.virx = 0.0;
406	–	virial.viry = 0.0;
407	–	virial.virz = 0.0;
245
246		for (i=1; i <= natom; i++)
247		atom[i].energy = 0.0;
#	Line 462 \| Line 299
299		{
300		double etot;
301		int i;
465	–	// struct timeb start,end;
302
303		energies.total = 0.0;
304		energies.estr = 0.0;
#	Line 549 \| Line 385
385		}
386
387		/* ================================================================== */
552	–
553	–	int iscoord_bond(int i, int j)
554	–	{
555	–	int k;
556	–	long int mask;
557	–
558	–	mask = 1 << METCOORD_MASK ;
559	–
560	–	for (k=0; k < MAXIAT; k++)
561	–	{
562	–	if (atom[i].iat[k] == j && atom[i].bo[k] == 9)
563	–	{
564	–	if (atom[i].type >= 300)
565	–	{
566	–	if (atom[j].type == 2 \|\| atom[j].type == 4 \|\| atom[j].type == 29 \|\|
567	–	atom[j].type == 30 \|\| atom[j].type == 40 \|\| atom[j].type == 48 )
568	–	return TRUE;
569	–	}
570	–	if (atom[j].type >= 300)
571	–	{
572	–	if (atom[i].type == 2 \|\| atom[i].type == 4 \|\| atom[i].type == 29 \|\|
573	–	atom[i].type == 30 \|\| atom[i].type == 40 \|\| atom[i].type == 48 )
574	–	return TRUE;
575	–	}
576	–	}
577	–	}
578	–	// throughout Metal-lone pairs if coordinated bond
579	–	// metal to donor atom
580	–	if ( (atom[i].type >= 300) && atom[j].flags & mask)
581	–	return TRUE;
582	–	if ( (atom[j].type >= 300) && atom[i].flags & mask)
583	–	return TRUE;
584	–
585	–	return FALSE;
586	–	}
587	–
588	–	void get_coordbonds()
589	–	{
590	–	int i, j, k, iatype, jatm;
591	–
592	–	cbonds.nbnd = 0;
593	–	pot.use_coordb = FALSE;
594	–
595	–	for (i=1; i <= natom; i++)
596	–	{
597	–	if (atom[i].mmx_type >= 300)
598	–	{
599	–	for (j=0; j <= MAXIAT; j++)
600	–	{
601	–	if (atom[i].bo[j] == 9) // coord bond
602	–	{
603	–	iatype = atom[atom[i].iat[j]].mmx_type;
604	–	if (iatype == 2 \|\| iatype == 4 \|\| iatype == 29 \|\|
605	–	iatype == 30 \|\| iatype == 40 \|\| iatype == 48 )
606	–	{
607	–	cbonds.icbonds[cbonds.nbnd][0] = i;
608	–	cbonds.icbonds[cbonds.nbnd][1] = atom[i].iat[j];
609	–	cbonds.icbonds[cbonds.nbnd][2] = 0;
610	–	cbonds.nbnd++;
611	–	} else
612	–	{
613	–	jatm = atom[i].iat[j];
614	–	for (k = 0; k < MAXIAT; k++)
615	–	{
616	–	if (atom[atom[jatm].iat[k]].mmx_type == 20) // lp
617	–	{
618	–	cbonds.icbonds[cbonds.nbnd][0] = i;
619	–	cbonds.icbonds[cbonds.nbnd][1] = jatm;
620	–	cbonds.icbonds[cbonds.nbnd][2] = atom[jatm].iat[k];
621	–	cbonds.nbnd++;
622	–	}
623	–	}
624	–	}
625	–	if (cbonds.nbnd > MAXCBND)
626	–	{
627	–	fprintf(pcmlogfile,"Error - Too many coordinated bonds!\nProgram will now quit.\n");
628	–	exit(0);
629	–	}
630	–	}
631	–	}
632	–	}
633	–	}
634	–	if (cbonds.nbnd > 0)
635	–	pot.use_coordb = TRUE;
636	–	}
388		int isbond(int i, int j)
389		{
390		int k;
#	Line 644 \| Line 395
395		}
396		return FALSE;
397		}
398	<
398	>	// ==========================
399		void get_bonds()
400		{
401		int i, j;
651	–	long int mask;
402
403	<	bonds_ff.nbnd = 0;
654	<	mask = 1L << 0; // flag for pi atoms
655	<
403	>	bonds_ff.nbnd = 0;
404		for (i=1; i <= natom; i++)
405		{
406		for (j=i+1; j <= natom; j++)

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines