PageRenderTime 61ms CodeModel.GetById 24ms RepoModel.GetById 0ms app.codeStats 0ms

/apbs-1.3-source/tools/manip/coulomb.c

#
C | 175 lines | 148 code | 18 blank | 9 comment | 23 complexity | 196473a53c493666ce422b9ea1312037 MD5 | raw file
Possible License(s): LGPL-2.0, BSD-3-Clause, GPL-2.0
  1. /**
  2. * @file coulomb.c
  3. * @author Nathan Baker
  4. * @brief Small program to calculate Coulombic energies
  5. * @version $Id: coulomb.c 1615 2010-10-20 19:16:35Z sobolevnrm $
  6. */
  7. #include "apbscfg.h"
  8. #include "maloc/maloc.h"
  9. #include "apbs/apbs.h"
  10. int main(int argc, char **argv) {
  11. /* OBJECTS */
  12. Valist *alist = VNULL;
  13. Vatom *atom = VNULL;
  14. Vgreen *green = VNULL;
  15. Vio *sock = VNULL;
  16. double *pos, energy, zmagic, disp[3], force[3];
  17. double *fx, *fy, *fz, *pot, *xp, *yp, *zp, *qp;
  18. int i, j;
  19. int doenergy = 0;
  20. int doforce = 0;
  21. /* SYSTEM PARAMETERS */
  22. char *path;
  23. char *usage = "\n\n\
  24. This program calculates electrostatic properties using Coulomb's law;\n\
  25. i.e., the Green's function for the Poisson operator in a homogeneous\n\
  26. medium with dielectric constant 1 (vacumm). It is very important to\n\
  27. realize that all energies, forces, etc. are calculated with this\n\
  28. dielectric constant of 1 and must be scaled to be compared with other\n\
  29. calculations.\n\n\
  30. Usage: coulomb [-e] [-f] <molecule.pqr>\n\n\
  31. where <molecule.pqr> is the path to the molecule\n\
  32. structure in PQR format. By default the total energies and forces\n\
  33. will be printed. This program also supports the following options:\n\
  34. -e give per-atom energies\n\
  35. -f give per-atom forces\n\n";
  36. Vio_start();
  37. if ((argc > 4) || (argc < 2)) {
  38. Vnm_print(2, "\n*** Syntax error: got %d arguments, expected 2.\n",
  39. argc);
  40. Vnm_print(2, "%s", usage);
  41. exit(666);
  42. };
  43. if (argc > 2) {
  44. for (i=1; i<(argc-1); i++) {
  45. if (strcmp("-e", argv[i]) == 0) {
  46. Vnm_print(1, "Providing per-atom energies...\n");
  47. doenergy = 1;
  48. } else if (strcmp("-f", argv[i]) == 0) {
  49. Vnm_print(1, "Providing per-atom forces...\n");
  50. doforce = 1;
  51. } else if (strcmp("-ef", argv[i]) == 0 || \
  52. strcmp("-fe", argv[i]) == 0){
  53. Vnm_print(1, "Providing per-atom forces and energies...\n");
  54. doforce = 1;
  55. doenergy = 1;
  56. } else {
  57. Vnm_print(2, "Ignoring option %s\n", argv[i]);
  58. }
  59. }
  60. path = argv[argc-1];
  61. } else {
  62. path = argv[1];
  63. }
  64. Vnm_print(1, "Setting up atom list from %s.\n", path);
  65. alist = Valist_ctor();
  66. sock = Vio_ctor("FILE", "ASC", VNULL, path, "r");
  67. if (sock == VNULL) {
  68. Vnm_print(2, "Problem opening virtual socket %s!\n",
  69. path);
  70. return 0;
  71. }
  72. if (Vio_accept(sock, 0) < 0) {
  73. Vnm_print(2, "Problem accepting virtual socket %s!\n",
  74. path);
  75. return 0;
  76. }
  77. Valist_readPQR(alist,VNULL,sock);
  78. Vnm_print(1, "Read %d atoms\n", Valist_getNumberAtoms(alist));
  79. Vnm_print(1, "Setting up Green's function object.\n");
  80. green = Vgreen_ctor(alist);
  81. /* Initialize variables */
  82. Vnm_print(1, "Dielectric constant = 1 (vacuum)\n");
  83. Vnm_print(1, "Distances in Angstroms\n");
  84. Vnm_print(1, "Charges in electrons\n");
  85. zmagic = (1e-3)*Vunit_ec*Vunit_Na;
  86. energy = 0.0;
  87. force[0] = 0.0; force[1] = 0.0; force[2] = 0.0;
  88. Vnm_print(1, "Allocating space for solution...\n");
  89. fx = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  90. fy = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  91. fz = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  92. pot = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  93. xp = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  94. yp = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  95. zp = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  96. qp = Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
  97. for (i=0; i<Valist_getNumberAtoms(alist); i++) {
  98. fx[i] = 0.0;
  99. fy[i] = 0.0;
  100. fz[i] = 0.0;
  101. pot[i] = 0.0;
  102. atom = Valist_getAtom(alist, i);
  103. pos = Vatom_getPosition(atom);
  104. xp[i] = pos[0];
  105. yp[i] = pos[1];
  106. zp[i] = pos[2];
  107. qp[i] = Vatom_getCharge(atom);
  108. }
  109. Vnm_print(1, "Calculating...\n");
  110. Vgreen_coulombD(green, Valist_getNumberAtoms(alist), xp, yp, zp,
  111. pot, fx, fy, fz);
  112. for (i=0; i<Valist_getNumberAtoms(alist); i++) {
  113. fx[i] *= (-0.5*qp[i]*zmagic);
  114. fy[i] *= (-0.5*qp[i]*zmagic);
  115. fz[i] *= (-0.5*qp[i]*zmagic);
  116. pot[i] *= (0.5*qp[i]*zmagic);
  117. energy += pot[i];
  118. force[0] += fx[i];
  119. force[1] += fy[i];
  120. force[2] += fz[i];
  121. if (doenergy) {
  122. Vnm_print(1, "\tAtom %d: Energy = %1.12E kJ/mol\n", i+1, pot[i]);
  123. }
  124. if (doforce) {
  125. Vnm_print(1, "\tAtom %d: x-force = %1.12E kJ/mol/A\n", i+1,
  126. fx[i]);
  127. Vnm_print(1, "\tAtom %d: y-force = %1.12E kJ/mol/A\n", i+1,
  128. fy[i]);
  129. Vnm_print(1, "\tAtom %d: z-force = %1.12E kJ/mol/A\n", i+1,
  130. fz[i]);
  131. }
  132. }
  133. Vnm_print(1, "\n\n-------------------------------------------------------\n");
  134. Vnm_print(1, "Total energy = %1.12e kJ/mol in vacuum.\n", energy);
  135. Vnm_print(1, "Total x-force = %1.12e kJ/mol/A in vacuum.\n", force[0]);
  136. Vnm_print(1, "Total y-force = %1.12e kJ/mol/A in vacuum.\n", force[1]);
  137. Vnm_print(1, "Total z-force = %1.12e kJ/mol/A in vacuum.\n", force[2]);
  138. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  139. (void **)&fx);
  140. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  141. (void **)&fy);
  142. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  143. (void **)&fz);
  144. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  145. (void **)&pot);
  146. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  147. (void **)&xp);
  148. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  149. (void **)&yp);
  150. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  151. (void **)&zp);
  152. Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double),
  153. (void **)&qp);
  154. Vgreen_dtor(&green);
  155. Valist_dtor(&alist);
  156. return 0;
  157. }