• R/O
  • SSH
  • HTTPS

molby: Commit


Commit MetaInfo

Revision69 (tree)
Time2011-07-06 02:46:20
Authortoshinagata1964

Log Message

Molecule#{bond,angle,dihedral,improper,vdw}_par were wrongly implemented. Fixed. 'Create SANDER input' is implemented (still not tested much).

Change Summary

Incremental Difference

--- trunk/Scripts/md.rb (revision 68)
+++ trunk/Scripts/md.rb (revision 69)
@@ -665,6 +665,355 @@
665665 }
666666 end
667667
668+ # Create sander input files from current molecule
669+ def export_prmtop
670+
671+ def error_dialog(msg)
672+ Dialog.run("AMBER Export Error") {
673+ layout(1, item(:text, :title=>msg))
674+ set_attr(1, :hidden=>true)
675+ }
676+ end
677+
678+ def format_print(fp, num, fmt, ary)
679+ fmts = "%#{fmt}%s"
680+ ary.each_with_index { |x, i|
681+ fp.printf(fmts, x, (i % num == num - 1 ? "\n" : ""))
682+ }
683+ fp.print "\n"
684+ end
685+
686+ par = self.parameter
687+
688+ # Create residue number table
689+ restable = []
690+ resno = -1
691+ self.each_atom { |ap|
692+ if ap.res_seq > resno
693+ restable.push(ap.index)
694+ resno = ap.res_seq
695+ elsif ap.res_seq < resno
696+ error_dialog("The residue number is not in ascending order. Please use 'sort by residue' command before creating AMBER input.")
697+ return nil
698+ end
699+ }
700+
701+ # Check if periodic box is used
702+ periodic = (self.box && self.box[4].all? { |n| n != 0 })
703+ if periodic
704+ fragments = []
705+ last_solute = nil
706+ first_solv_mol = nil
707+ self.each_fragment { |gr|
708+ if gr.range_at(1) != nil
709+ msg = gr.to_s.sub(/IntGroup/, "")
710+ error_dialog("The atoms #{msg} are one fragment, but the atom indices are not consecutive.")
711+ return nil
712+ end
713+ fragments.push(gr.length)
714+ if last_solute == nil && self.atoms[gr[0]].seq_name == "SOLV"
715+ # Calculate the residue number of the last solute atom (gr[0] - 1)
716+ n = gr[0] - 1
717+ if n < 0
718+ error_dialog("All atoms are labeled as the solvent??")
719+ return nil
720+ end
721+ first_solv_mol = fragments.length
722+ restable.each_with_index { |m, i|
723+ n -= m
724+ if n <= 0
725+ last_solute = i + 1
726+ break
727+ end
728+ }
729+ end
730+ }
731+ if last_solute == nil
732+ last_solute = restable.length
733+ first_solv_mol = fragments.length
734+ end
735+ end
736+
737+ # Count bonds and angles, including and not-including hydrogens
738+ bonds_h = []
739+ bonds_a = []
740+ self.bonds.each_with_index { |e, i|
741+ if e.any? { |n| self.atoms[n].atomic_number == 1 }
742+ bonds_h.push(i)
743+ else
744+ bonds_a.push(i)
745+ end
746+ }
747+ angles_h = []
748+ angles_a = []
749+ self.angles.each_with_index { |e, i|
750+ if e.any? { |n| self.atoms[n].atomic_number == 1 }
751+ angles_h.push(i)
752+ else
753+ angles_a.push(i)
754+ end
755+ }
756+
757+ # Build exclusion table (before processing dihedral angles)
758+ exnumbers = []
759+ extable = []
760+ exhash = {}
761+ self.each_atom { |ap|
762+ ex = ap.exclusion
763+ exx = []
764+ ex.each_with_index { |x, i|
765+ x.each { |n|
766+ if n > ap.index
767+ exhash[ap.index * 1000000 + n] = i + 2 # 2 for 1-2, 3 for 1-3, 4 for 1-4
768+ exx.push(n + 1)
769+ end
770+ }
771+ }
772+ if exx.length > 0
773+ exx.sort!
774+ else
775+ exx.push(0)
776+ end
777+ exnumbers.push(exx.length)
778+ extable.concat(exx)
779+ }
780+
781+ # Count dihedrals and impropers, including and not-including hydrogens
782+ dihedrals_h = []
783+ dihedrals_a = []
784+ self.dihedrals.each_with_index { |e, i|
785+ flag = 0
786+ k = (e[0] < e[3] ? e[0] * 1000000 + e[3] : e[3] * 1000000 + e[0])
787+ if exhash[k] == 4
788+ exhash[k] = -4 # Handle 1,4-exclusion for only one dihedrals
789+ elsif exhash[k] == -4
790+ flag = 1000000 # Skip calculation of 1,4-interaction to avoid double counting
791+ elsif exhash[k] && exhash[k] > 0
792+ flag = 1000000 # 5-member or smaller ring: e[0]...e[3] is 1-2 or 1-3 exclusion pair
793+ end
794+ if e.any? { |n| self.atoms[n].atomic_number == 1 }
795+ dihedrals_h.push(i + flag)
796+ else
797+ dihedrals_a.push(i + flag)
798+ end
799+ }
800+ self.impropers.each_with_index { |e, i|
801+ if e.any? { |n| self.atoms[n].atomic_number == 1 }
802+ dihedrals_h.push(i + 2000000)
803+ else
804+ dihedrals_a.push(i + 2000000)
805+ end
806+ }
807+
808+ # Create vdw parameter table
809+ nvdw_pars = par.nvdws
810+ vdw_pair_table = []
811+ vdw_access_table = []
812+ (0..nvdw_pars - 1).each { |i|
813+ p1 = par.vdws[i]
814+ (0..i).each { |j|
815+ p2 = par.vdws[j]
816+ pp = par.lookup("vdw_pair", [p1.atom_type, p2.atom_type])
817+ if pp
818+ # Vdw pair parameter is defined
819+ eps = pp.eps
820+ d = pp.r_eq * 2
821+ else
822+ eps = Math.sqrt(p1.eps * p2.eps)
823+ d = p1.r_eq + p2.r_eq
824+ end
825+ vdw_access_table[i * nvdw_pars + j] = vdw_access_table[j * nvdw_pars + i] = vdw_pair_table.length
826+ vdw_pair_table.push([(d ** 12) * eps, 2 * (d ** 6) * eps])
827+ }
828+ }
829+
830+ basename = (self.path ? File.basename(self.path, ".*") : self.name)
831+ fname = Dialog.save_panel("AMBER prmtop/inpcrd file name", self.dir, basename + ".prmtop", "All files|*.*")
832+ return nil if !fname
833+
834+ open(fname, "w") { |fp|
835+ date = Time.now.localtime.strftime("%m/%d/%y %H:%M:%S")
836+
837+ fp.print "%VERSION VERSION_STAMP = V0001.000 DATE = #{date}\n"
838+
839+ fp.print "%FLAG TITLE\n%FORMAT(20a4)\n"
840+ fp.printf "%-80.80s\n", self.name
841+
842+ fp.print "%FLAG POINTERS\n%FORMAT(10I8)\n"
843+ fp.printf "%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", self.natoms, par.nvdws, bonds_h.length, bonds_a.length, angles_h.length, angles_a.length, dihedrals_h.length, dihedrals_a.length, 0, 0
844+ fp.printf "%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", extable.length, restable.length, bonds_a.length, angles_a.length, dihedrals_a.length, par.nbonds, par.nangles, par.ndihedrals + par.nimpropers, par.nvdws, 0
845+ fp.printf "%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n", 0, 0, 0, 0, 0, 0, 0, (periodic ? 1 : 0), restable.max, 0
846+ fp.printf "%8d\n", 0
847+
848+ fp.print "%FLAG ATOM_NAME\n%FORMAT(20a4)\n"
849+ format_print(fp, 20, "-4.4s", self.atoms.map { |ap| ap.name })
850+
851+ fp.print "%FLAG CHARGE\n%FORMAT(5E16.8)\n"
852+ format_print(fp, 5, "16.8E", self.atoms.map { |ap| ap.charge * 18.2223 })
853+
854+ fp.print "%FLAG MASS\n%FORMAT(5E16.8)\n"
855+ format_print(fp, 5, "16.8E", self.atoms.map { |ap| ap.weight })
856+
857+ fp.print "%FLAG ATOM_TYPE_INDEX\n%FORMAT(10I8)\n"
858+ format_print(fp, 10, "8d", (0...self.natoms).map { |i| self.vdw_par(i).index + 1 })
859+
860+ fp.print "%FLAG NUMBER_EXCLUDED_ATOMS\n%FORMAT(10I8)\n"
861+ format_print(fp, 10, "8d", exnumbers)
862+
863+ fp.print "%FLAG NONBONDED_PARM_INDEX\n%FORMAT(10I8)\n"
864+ format_print(fp, 10, "8d", vdw_access_table.map { |n| n + 1 })
865+
866+ fp.print "%FLAG RESIDUE_LABEL\n%FORMAT(20a4)\n"
867+ format_print(fp, 20, "-4.4s", restable.map { |n| self.atoms[n].res_name })
868+
869+ fp.print "%FLAG RESIDUE_POINTER\n%FORMAT(10I8)\n"
870+ format_print(fp, 10, "8d", restable.map { |n| n + 1 })
871+
872+ fp.print "%FLAG BOND_FORCE_CONSTANT\n%FORMAT(5E16.8)\n"
873+ format_print(fp, 5, "16.8E", par.bonds.map { |p| p.k })
874+
875+ fp.print "%FLAG BOND_EQUIL_VALUE\n%FORMAT(5E16.8)\n"
876+ format_print(fp, 5, "16.8E", par.bonds.map { |p| p.r0 })
877+
878+ fp.print "%FLAG ANGLE_FORCE_CONSTANT\n%FORMAT(5E16.8)\n"
879+ format_print(fp, 5, "16.8E", par.angles.map { |p| p.k })
880+
881+ fp.print "%FLAG ANGLE_EQUIL_VALUE\n%FORMAT(5E16.8)\n"
882+ format_print(fp, 5, "16.8E", par.angles.map { |p| p.a0 * Math::PI / 180.0 })
883+
884+ fp.print "%FLAG DIHEDRAL_FORCE_CONSTANT\n%FORMAT(5E16.8)\n"
885+ format_print(fp, 5, "16.8E", par.dihedrals.map { |p| p.k } + par.impropers.map { |p| p.k })
886+
887+ fp.print "%FLAG DIHEDRAL_PERIODICITY\n%FORMAT(5E16.8)\n"
888+ format_print(fp, 5, "16.8E", par.dihedrals.map { |p| p.period } + par.impropers.map { |p| p.period })
889+
890+ fp.print "%FLAG DIHEDRAL_PHASE\n%FORMAT(5E16.8)\n"
891+ format_print(fp, 5, "16.8E", par.dihedrals.map { |p| p.phi0 * Math::PI / 180.0 } + par.impropers.map { |p| p.phi0 * Math::PI / 180.0 })
892+
893+ fp.print "%FLAG SCEE_SCALE_FACTOR\n%FORMAT(5E16.8)\n"
894+ format_print(fp, 5, "16.8E", par.dihedrals.map { |p| 1.2 } + par.impropers.map { |p| 0.0 })
895+
896+ fp.print "%FLAG SCNB_SCALE_FACTOR\n%FORMAT(5E16.8)\n"
897+ format_print(fp, 5, "16.8E", par.dihedrals.map { |p| 2.0 } + par.impropers.map { |p| 0.0 })
898+
899+ fp.print "%FLAG SOLTY\n%FORMAT(5E16.8)\n"
900+ format_print(fp, 5, "16.8E", (0...par.nvdws).map { 0.0 } )
901+
902+ fp.print "%FLAG LENNARD_JONES_ACOEF\n%FORMAT(5E16.8)\n"
903+ format_print(fp, 5, "16.8E", vdw_pair_table.map { |x| x[0] })
904+
905+ fp.print "%FLAG LENNARD_JONES_BCOEF\n%FORMAT(5E16.8)\n"
906+ format_print(fp, 5, "16.8E", vdw_pair_table.map { |x| x[1] })
907+
908+ fp.print "%FLAG BONDS_INC_HYDROGEN\n%FORMAT(10I8)\n"
909+ format_print(fp, 10, "8d", bonds_h.map { |n|
910+ x = self.bonds[n]
911+ [x[0] * 3, x[1] * 3, self.bond_par(n).index + 1] }.flatten)
912+
913+ fp.print "%FLAG BONDS_WITHOUT_HYDROGEN\n%FORMAT(10I8)\n"
914+ format_print(fp, 10, "8d", bonds_a.map { |n|
915+ x = self.bonds[n]
916+ [x[0] * 3, x[1] * 3, self.bond_par(n).index + 1] }.flatten)
917+
918+ fp.print "%FLAG ANGLES_INC_HYDROGEN\n%FORMAT(10I8)\n"
919+ format_print(fp, 10, "8d", angles_h.map { |n|
920+ x = self.angles[n]
921+ [x[0] * 3, x[1] * 3, x[2] * 3, self.angle_par(n).index + 1] }.flatten)
922+
923+ fp.print "%FLAG ANGLES_WITHOUT_HYDROGEN\n%FORMAT(10I8)\n"
924+ format_print(fp, 10, "8d", angles_a.map { |n|
925+ x = self.angles[n]
926+ [x[0] * 3, x[1] * 3, x[2] * 3, self.angle_par(n).index + 1] }.flatten)
927+
928+ [dihedrals_h, dihedrals_a].each { |dihed|
929+ if dihed == dihedrals_h
930+ fp.print "%FLAG DIHEDRALS_INC_HYDROGEN\n"
931+ else
932+ fp.print "%FLAG DIHEDRALS_WITHOUT_HYDROGEN\n"
933+ end
934+ fp.print "%FORMAT(10I8)\n"
935+ format_print(fp, 10, "8d", dihed.map { |n|
936+ if n < 2000000
937+ x = self.dihedrals[n % 1000000]
938+ # Note: if n >= 1000000, then the 1-4 interaction should be ignored to avoid double calculation. This is specified by the negative index for the third atom, but if the third atom is 0 we have a problem. To avoid this situation, the atom array is reversed.
939+ if n >= 1000000 && x[2] == 0
940+ x = x.reverse
941+ end
942+ k = self.dihedral_par(n % 1000000).index + 1
943+ [x[0] * 3, x[1] * 3, (n >= 1000000 ? -x[2] : x[2]) * 3, x[3] * 3, k]
944+ else
945+ x = self.impropers[n % 1000000]
946+ # The improper torsion is specified by the negative index for the fourth atom, but if the fourth atom is 0 then we have a problem. To avoid this situation, the first and fourth atoms are exchanged.
947+ if x[3] == 0
948+ x = [x[3], x[1], x[2], x[0]]
949+ end
950+ k = self.improper_par(n % 1000000).index + 1 + par.ndihedrals
951+ [x[0] * 3, x[1] * 3, -x[2] * 3, -x[3] * 3, k]
952+ end
953+ }.flatten)
954+ } # end each [dihedral_h, dihedral_a]
955+
956+ fp.print "%FLAG EXCLUDED_ATOMS_LIST\n%FORMAT(10I8)\n"
957+ format_print(fp, 10, "8d", extable)
958+
959+ fp.print "%FLAG HBOND_ACOEF\n%FORMAT(5E16.8)\n"
960+ fp.print "\n"
961+
962+ fp.print "%FLAG HBOND_BCOEF\n%FORMAT(5E16.8)\n"
963+ fp.print "\n"
964+
965+ fp.print "%FLAG HBCUT\n%FORMAT(5E16.8)\n"
966+ fp.print "\n"
967+
968+ fp.print "%FLAG AMBER_ATOM_TYPE\n%FORMAT(20a4)\n"
969+ format_print(fp, 20, "-4.4s", self.atoms.map { |ap| ap.atom_type })
970+
971+ fp.print "%FLAG TREE_CHAIN_CLASSIFICATION\n%FORMAT(20a4)\n"
972+ format_print(fp, 20, "-4.4s", self.atoms.map { |ap| "M" })
973+
974+ fp.print "%FLAG JOIN_ARRAY\n%FORMAT(10I8)\n"
975+ format_print(fp, 10, "8d", self.atoms.map { |ap| 0 })
976+
977+ fp.print "%FLAG IROTAT\n%FORMAT(10I8)\n"
978+ format_print(fp, 10, "8d", self.atoms.map { |ap| 0 })
979+
980+ fp.print "%FLAG RADIUS_SET\n%FORMAT(1a80)\n"
981+ fp.print "modified Bondi radii (mbondi)\n"
982+
983+ fp.print "%FLAG RADII\n%FORMAT(5E16.8)\n"
984+ format_print(fp, 5, "16.8E", self.atoms.map { |ap|
985+ (ap.atomic_number == 1 ? 1.30 : 1.70) })
986+
987+ fp.print "%FLAG SCREEN\n%FORMAT(5E16.8)\n"
988+ format_print(fp, 5, "16.8E", self.atoms.map { |ap|
989+ (ap.atomic_number == 1 ? 0.85 : 0.72) })
990+
991+ if periodic
992+ fp.print "%FLAG SOLVENT_POINTERS\n%FORMAT(3I8)\n"
993+ fp.printf "%8d%8d%8d\n", last_solute, fragments.length, first_solv_mol
994+
995+ fp.print "%FLAG ATOMS_PER_MOLECULE\n%FORMAT(10I8)\n"
996+ format_print(fp, 10, "8d", fragments)
997+
998+ fp.print "%FLAG BOX_DIMENSIONS\n%FORMAT(5E16.8)\n"
999+ box = self.box
1000+ fp.printf "%16.8E%16.8E%16.8E%16.8E\n", 0, box[0].length, box[1].length, box[2].length
1001+ end
1002+ }
1003+
1004+ fname2 = fname.sub(/\.prmtop$/, '.inpcrd')
1005+ fname2 += '.inpcrd' if fname == fname2
1006+
1007+ open(fname2, "w") { |fp|
1008+ fp.printf "%-80.80s\n", self.name
1009+ fp.printf "%6d\n", self.natoms
1010+ format_print(fp, 6, "12.7f", self.atoms.map { |ap| [ap.x, ap.y, ap.z] }.flatten )
1011+ }
1012+
1013+ return true
1014+
1015+ end
1016+
6681017 def count_elements
6691018 elements = []
6701019 each_atom { |ap|
--- trunk/wxSources/MyDocument.cpp (revision 68)
+++ trunk/wxSources/MyDocument.cpp (revision 69)
@@ -99,6 +99,7 @@
9999 EVT_MENU(myMenuID_ExpandBySymmetry, MyDocument::OnExpandBySymmetry)
100100 EVT_MENU(myMenuID_RunAntechamber, MyDocument::OnInvokeAntechamber)
101101 EVT_MENU(myMenuID_RunResp, MyDocument::OnInvokeResp)
102+ EVT_MENU(myMenuID_CreateSanderInput, MyDocument::OnCreateSanderInput)
102103 EVT_MENU(myMenuID_CreateGamessInput, MyDocument::OnCreateGamessInput)
103104 EVT_MENU(myMenuID_CreateMOCube, MyDocument::OnCreateMOCube)
104105 EVT_MENU(myMenuID_ShowAllAtoms, MyDocument::OnShowAllAtoms)
@@ -869,6 +870,49 @@
869870 {
870871 }
871872
873+static wxString
874+sCreateTemporaryLogDirectoryForAC(const wxString& filename)
875+{
876+// char *ante_dir;
877+ char *log_dir;
878+ int i, status;
879+
880+ /* Extract the name */
881+ wxFileName fname(filename);
882+ wxString name = fname.GetName();
883+
884+ status = MyAppCallback_getGlobalSettingsWithType("antechamber.log_dir", 's', &log_dir);
885+ if (status) {
886+ char *hdir = MyAppCallback_getDocumentHomeDir();
887+ asprintf(&log_dir, "%s/antechamber", (hdir ? hdir : ""));
888+ if (hdir != NULL)
889+ free(hdir);
890+ }
891+ fix_dosish_path(log_dir);
892+
893+ wxString tdir;
894+
895+ /* Prepare the log directory */
896+ wxString dirname(log_dir, wxConvFile);
897+ if (!wxFileName::Mkdir(dirname, 0777, wxPATH_MKDIR_FULL)) {
898+ MyAppCallback_errorMessageBox("Cannot create log directory '%s'", log_dir);
899+ free(log_dir);
900+ return tdir; /* empty */
901+ }
902+ free(log_dir);
903+
904+ for (i = 0; i < 1000; i++) {
905+ tdir = dirname + wxFileName::GetPathSeparator() + name + wxString::Format(_T("_%04d"), i);
906+ if (!wxFileName::DirExists(tdir))
907+ break;
908+ }
909+ if (i >= 1000 || !wxFileName::Mkdir(tdir)) {
910+ MyAppCallback_errorMessageBox("Cannot create temporary files. Please make sure the log directory has enough space for writing.");
911+ tdir.Empty();
912+ }
913+ return tdir;
914+}
915+
872916 static bool
873917 sRemoveDirectoryRecursively(const wxString &dir)
874918 {
@@ -900,53 +944,111 @@
900944 return ::wxRmdir(dir);
901945 }
902946
947+static int
948+sEraseLogFiles(const wxString& tdir, int status)
949+{
950+ bool success = true;
951+ Int log_keep_number, n, i, j;
952+ char *log_level;
953+ wxString dir2;
954+
955+ if (MyAppCallback_getGlobalSettingsWithType("antechamber.log_level", 's', &log_level) != 0)
956+ log_level = NULL;
957+ if (MyAppCallback_getGlobalSettingsWithType("antechamber.log_keep_number", 'i', &log_keep_number) != 0)
958+ log_keep_number = 5;
959+ if (log_level == NULL || strcmp(log_level, "none") == 0 || (strcmp(log_level, "error_only") == 0 && status == 0)) {
960+ // Erase the present log
961+ if (!sRemoveDirectoryRecursively(tdir)) {
962+ success = false;
963+ dir2 = tdir;
964+ }
965+ } else if (strcmp(log_level, "latest") == 0) {
966+ wxString dirname = tdir.BeforeLast(wxFileName::GetPathSeparator());
967+ wxDir wdir(dirname);
968+ wxString name;
969+ wxArrayString files;
970+ n = 0;
971+ if (wdir.GetFirst(&name)) {
972+ do {
973+ wxString fullname = dirname + wxFileName::GetPathSeparator() + name;
974+ if (wxDir::Exists(fullname)) {
975+ files.Add(fullname);
976+ n++;
977+ }
978+ } while (wdir.GetNext(&name));
979+ }
980+ if (n > log_keep_number) {
981+ // Sort directories by creation date
982+ struct temp_struct { time_t tm; int idx; } *tp;
983+ tp = (struct temp_struct *)malloc(sizeof(struct temp_struct) * n);
984+ for (i = 0; i < n; i++) {
985+ wxFileName fn(files[i], wxEmptyString);
986+ wxDateTime dt;
987+ j = fn.GetTimes(NULL, NULL, &dt);
988+ tp[i].tm = dt.GetTicks();
989+ tp[i].idx = i;
990+ }
991+ for (i = 0; i < n; i++) {
992+ struct temp_struct temp;
993+ int k = i;
994+ for (j = i + 1; j < n; j++) {
995+ if (tp[j].tm < tp[k].tm)
996+ k = j;
997+ }
998+ if (k != i) {
999+ temp = tp[k];
1000+ tp[k] = tp[i];
1001+ tp[i] = temp;
1002+ }
1003+ }
1004+ // Keep last log_keep_number and delete the rest
1005+ for (i = 0; i < n - log_keep_number; i++) {
1006+ if (!sRemoveDirectoryRecursively(files[tp[i].idx])) {
1007+ success = false;
1008+ dir2 = files[tp[i].idx];
1009+ break;
1010+ }
1011+ }
1012+ }
1013+ }
1014+
1015+ if (success) {
1016+ return 0;
1017+ } else {
1018+ MyAppCallback_errorMessageBox("Error during deleting log file '%s'", (const char *)dir2.mb_str(wxConvUTF8));
1019+ return -1;
1020+ }
1021+}
1022+
9031023 void
9041024 MyDocument::OnInvokeAntechamber(wxCommandEvent &event)
9051025 {
906- char *ante_dir;
907- char *log_dir, *log_level, buf[256];
908- Int net_charge, i, j, n, log_keep_number, calc_charge, use_residue;
1026+ char *ante_dir, buf[256];
1027+ Int net_charge, i, n, calc_charge, use_residue;
9091028 int status;
9101029
9111030 /* Find the ambertool directory */
9121031 wxString ante = MyApp::FindResourcePath() + wxFileName::GetPathSeparator() + _T("amber11") + wxFileName::GetPathSeparator() + _T("bin");
9131032 ante_dir = strdup(ante.mb_str(wxConvFile));
1033+ fix_dosish_path(ante_dir);
9141034
9151035 /* Ask for antechamber options and log directory */
9161036 MolActionCreateAndPerform(mol, SCRIPT_ACTION(";i"), "cmd_antechamber", &n);
9171037 if (!n)
9181038 return;
1039+
9191040 if ((status = MyAppCallback_getGlobalSettingsWithType("antechamber.nc", 'i', &net_charge))
920- || (status = MyAppCallback_getGlobalSettingsWithType("antechamber.log_level", 's', &log_level))
921- || (status = MyAppCallback_getGlobalSettingsWithType("antechamber.log_keep_number", 'i', &log_keep_number))
9221041 || (status = MyAppCallback_getGlobalSettingsWithType("antechamber.calc_charge", 'i', &calc_charge))
923- || (status = MyAppCallback_getGlobalSettingsWithType("antechamber.use_residue", 'i', &use_residue))
924- || (status = MyAppCallback_getGlobalSettingsWithType("antechamber.log_dir", 's', &log_dir))) {
1042+ || (status = MyAppCallback_getGlobalSettingsWithType("antechamber.use_residue", 'i', &use_residue))) {
9251043 Molby_showError(status);
9261044 return;
9271045 }
928- fix_dosish_path(ante_dir);
929- fix_dosish_path(log_dir);
9301046
9311047 /* Prepare the log directory */
932- wxString dirname(log_dir, wxConvFile);
933- if (!wxFileName::Mkdir(dirname, 0777, wxPATH_MKDIR_FULL)) {
934- MyAppCallback_errorMessageBox("Cannot create log directory '%s'", log_dir);
1048+ wxString tdir = sCreateTemporaryLogDirectoryForAC(GetFilename());
1049+ if (tdir.IsEmpty())
9351050 return;
936- }
937- wxFileName filename(GetFilename());
938- wxString name = filename.GetName();
939- wxString tdir;
940- for (i = 0; i < 1000; i++) {
941- tdir = dirname + wxFileName::GetPathSeparator() + name + wxString::Format(_T("_%04d"), i);
942- if (!wxFileName::DirExists(tdir))
943- break;
944- }
945- if (i >= 1000 || !wxFileName::Mkdir(tdir)) {
946- MyAppCallback_errorMessageBox("Cannot create temporary files. Please make sure the log directory has enough space for writing.");
947- return;
948- }
949-
1051+
9501052 /* Move to the temporary directory and export the molecule as a pdb */
9511053 wxString cwd = wxFileName::GetCwd();
9521054 if (!wxFileName::SetCwd(tdir)) {
@@ -1048,64 +1150,7 @@
10481150 wxFileName::SetCwd(cwd);
10491151
10501152 /* Erase log files */
1051- bool success = true;
1052- wxString dir2;
1053- if (strcmp(log_level, "none") == 0 || (strcmp(log_level, "error_only") == 0 && status == 0)) {
1054- /* Erase the present log */
1055- if (!sRemoveDirectoryRecursively(tdir)) {
1056- success = false;
1057- dir2 = tdir;
1058- }
1059- } else if (strcmp(log_level, "latest") == 0) {
1060- wxDir wdir(dirname);
1061- wxString name;
1062- wxArrayString files;
1063- n = 0;
1064- if (wdir.GetFirst(&name)) {
1065- do {
1066- wxString fullname = dirname + wxFileName::GetPathSeparator() + name;
1067- if (wxDir::Exists(fullname)) {
1068- files.Add(fullname);
1069- n++;
1070- }
1071- } while (wdir.GetNext(&name));
1072- }
1073- if (n > log_keep_number) {
1074- /* Sort directories by creation date */
1075- struct temp_struct { time_t tm; int idx; } *tp;
1076- tp = (struct temp_struct *)malloc(sizeof(struct temp_struct) * n);
1077- for (i = 0; i < n; i++) {
1078- wxFileName fn(files[i], wxEmptyString);
1079- wxDateTime dt;
1080- j = fn.GetTimes(NULL, NULL, &dt);
1081- tp[i].tm = dt.GetTicks();
1082- tp[i].idx = i;
1083- }
1084- for (i = 0; i < n; i++) {
1085- struct temp_struct temp;
1086- int k = i;
1087- for (j = i + 1; j < n; j++) {
1088- if (tp[j].tm < tp[k].tm)
1089- k = j;
1090- }
1091- if (k != i) {
1092- temp = tp[k];
1093- tp[k] = tp[i];
1094- tp[i] = temp;
1095- }
1096- }
1097- /* Keep last log_keep_number and delete the rest */
1098- for (i = 0; i < n - log_keep_number; i++) {
1099- if (!sRemoveDirectoryRecursively(files[tp[i].idx])) {
1100- success = false;
1101- dir2 = files[tp[i].idx];
1102- break;
1103- }
1104- }
1105- }
1106- }
1107- if (!success)
1108- MyAppCallback_errorMessageBox("Error during deleting log file '%s'", (const char *)dir2.mb_str(wxConvUTF8));
1153+ sEraseLogFiles(tdir, status);
11091154
11101155 if (status == 0) {
11111156 ((MoleculeView *)GetFirstView())->GetListCtrl()->Update();
@@ -1120,6 +1165,12 @@
11201165 }
11211166
11221167 void
1168+MyDocument::OnCreateSanderInput(wxCommandEvent &event)
1169+{
1170+ MolActionCreateAndPerform(mol, SCRIPT_ACTION(""), "export_prmtop");
1171+}
1172+
1173+void
11231174 MyDocument::OnCreateGamessInput(wxCommandEvent &event)
11241175 {
11251176 MolActionCreateAndPerform(mol, SCRIPT_ACTION(""), "cmd_create_gamess_input");
--- trunk/wxSources/MyApp.h (revision 68)
+++ trunk/wxSources/MyApp.h (revision 69)
@@ -83,6 +83,7 @@
8383 myMenuID_MDTools = 151,
8484 myMenuID_RunAntechamber = 152,
8585 myMenuID_RunResp = 153,
86+ myMenuID_CreateSanderInput = 154,
8687 myMenuID_CreateGamessInput = 160,
8788 myMenuID_CreateMOCube = 161,
8889 myMenuID_ExecuteScript = 200,
--- trunk/wxSources/MyApp.cpp (revision 68)
+++ trunk/wxSources/MyApp.cpp (revision 69)
@@ -395,6 +395,7 @@
395395 wxMenu *md_tools_menu = new wxMenu;
396396 md_tools_menu->Append(myMenuID_RunAntechamber, _T("Antechamber/parmchk..."));
397397 md_tools_menu->Append(myMenuID_RunResp, _T("GAMESS/RESP..."));
398+ md_tools_menu->Append(myMenuID_CreateSanderInput, _T("Create SANDER input..."));
398399 md_menu->Append(myMenuID_MDTools, _T("Tools"), md_tools_menu);
399400
400401 wxMenu *qc_menu = new wxMenu;
@@ -1052,10 +1053,13 @@
10521053 int
10531054 MyAppCallback_getGlobalSettingsWithType(const char *key, int type, void *ptr)
10541055 {
1055- const char *s = MyAppCallback_getGlobalSettings(key);
1056+ int retval;
1057+ char *s = MyAppCallback_getGlobalSettings(key);
10561058 char desc[] = SCRIPT_ACTION("s; ");
10571059 desc[sizeof(desc) - 2] = type;
1058- return MolActionCreateAndPerform(NULL, desc, "eval", s, ptr);
1060+ retval = MolActionCreateAndPerform(NULL, desc, "eval", s, ptr);
1061+ free(s);
1062+ return retval;
10591063 }
10601064
10611065 int
--- trunk/wxSources/MyDocument.h (revision 68)
+++ trunk/wxSources/MyDocument.h (revision 69)
@@ -114,7 +114,8 @@
114114
115115 void OnInvokeResp(wxCommandEvent &event);
116116 void OnInvokeAntechamber(wxCommandEvent &event);
117-
117+ void OnCreateSanderInput(wxCommandEvent &event);
118+
118119 void OnCreateGamessInput(wxCommandEvent &event);
119120 void OnCreateMOCube(wxCommandEvent &event);
120121
--- trunk/memo.txt (revision 68)
+++ trunk/memo.txt (revision 69)
@@ -136,3 +136,42 @@
136136 The native format (mbsf) now preserves the display conditions, such as scale, orientation etc.
137137 Pasting MD parameters to a molecule with no parameters now works correctly.
138138 Importing GaussianW fch(k) files was not working; fixed. The extension ".fch" is now recognized as a Gaussian formatted checkpoint files.
139+
140+2011.7.5.
141+sander 用入力を作れるようにする。prmtop, inpcrd は一応作れるようになった。あとは sander の入力ファイル。関係ありそうな設定パラメータ:
142+
143+imin = 0 (no minimization), 1 (minimization)
144+ntx = 1 (read coords), 5 (read coords and velocities and box [if ntb>0]. velocities are used only if irest > 0.
145+irest = 0 (no restart), 1 (restart)
146+ntpr : every NTPR steps energy info are printed to mdout/mdinfo (default 50).
147+ntwx : every ntwx steps coords are written to mdcrd (default 0: no output)
148+ntf: force evaluation
149+ntb: 0 (no periodicity), 1 (constant volume), 2 (constant pressure)
150+cut: nonbonded cutoff, default 8.0 (for PME this is good); for igb > 0, larger values should be used
151+igb: 0 (no GB), 1 (Hawkins/Cramer/Truhlar GB), 2 (Onufriev/Bashford/Case GB), etc.
152+saltcon: salt concentraction (M), GB plus Debye-Huckel model.
153+rgbmax: cutoff for estimation of effective Born radii. Default 25 A.
154+gbsa: 0 (no surface area calc), 1 (LCPO surface area model), 2 (recursive approximation of spheres; used only in single point calc)
155+surften: surface tension. default 0.005 kcal/mol/A^2
156+rdt: (only used in LES)
157+
158+maxcyc: maximum number of cycles of minimization
159+nstlim: number of MD-steps to be performed
160+dt: time step (psec). Recommended value is 0.002 when SHAKE is used or 0.001 otherwise.
161+nrespa: frequency for calculation of slow forces. nrespa * dt > 0.004 causes instablitity.
162+
163+ntt: temperature control 0: constant energy, 1: constant temperature, 2: Andersen temperature, 3: Langevin dynamics
164+temp0: target temperature
165+tempi: initial temperature
166+tautp: time constant for heat bath coupling. Should be 0.5-5.0ps.
167+gamma_ln: collision frequency for Langevin dynamics.
168+vlimit: limiting velocity, default 20.
169+
170+ntp: constant pressure dynamics 0: no pressure scaling, 1: isotropic position scaling, 2: anisotropic pressure scaling
171+pres0: reference pressure (default 1 bar)
172+taup: pressure relaxation time (default 1.0 ps)
173+
174+ntc: SHAKE algorithm. generally ntf and ntc should be the same. 1: SHAKE is not perfomred, 2: bonds involving hydrogen are constrained, 3: all bonds are constrained
175+
176+
177+
--- trunk/MolLib/Ruby_bind/ruby_bind.c (revision 68)
+++ trunk/MolLib/Ruby_bind/ruby_bind.c (revision 69)
@@ -4694,6 +4694,8 @@
46944694 s_Molecule_BondPar(VALUE self, VALUE val)
46954695 {
46964696 Molecule *mol;
4697+ BondPar *bp;
4698+ UInt t1, t2;
46974699 Int ival;
46984700 Data_Get_Struct(self, Molecule, mol);
46994701 ival = NUM2INT(rb_Integer(val));
@@ -4702,7 +4704,12 @@
47024704 if (ival < 0)
47034705 ival += mol->nbonds;
47044706 s_RebuildMDParameterIfNecessary(self);
4705- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kBondParType, mol->arena->bond_par_i[ival]);
4707+ t1 = ATOM_AT_INDEX(mol->atoms, mol->bonds[ival * 2])->type;
4708+ t2 = ATOM_AT_INDEX(mol->atoms, mol->bonds[ival * 2 + 1])->type;
4709+ bp = ParameterLookupBondPar(mol->par, t1, t2, 0);
4710+ if (bp == NULL)
4711+ return Qnil;
4712+ return ValueFromMoleculeWithParameterTypeAndIndex(mol, kBondParType, bp - mol->par->bondPars);
47064713 }
47074714
47084715 /*
@@ -4715,6 +4722,8 @@
47154722 s_Molecule_AnglePar(VALUE self, VALUE val)
47164723 {
47174724 Molecule *mol;
4725+ AnglePar *ap;
4726+ UInt t1, t2, t3;
47184727 Int ival;
47194728 Data_Get_Struct(self, Molecule, mol);
47204729 ival = NUM2INT(rb_Integer(val));
@@ -4723,7 +4732,13 @@
47234732 if (ival < 0)
47244733 ival += mol->nangles;
47254734 s_RebuildMDParameterIfNecessary(self);
4726- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kAngleParType, mol->arena->angle_par_i[ival]);
4735+ t1 = ATOM_AT_INDEX(mol->atoms, mol->angles[ival * 3])->type;
4736+ t2 = ATOM_AT_INDEX(mol->atoms, mol->angles[ival * 3 + 1])->type;
4737+ t3 = ATOM_AT_INDEX(mol->atoms, mol->angles[ival * 3 + 2])->type;
4738+ ap = ParameterLookupAnglePar(mol->par, t1, t2, t3, 0);
4739+ if (ap == NULL)
4740+ return Qnil;
4741+ return ValueFromMoleculeWithParameterTypeAndIndex(mol, kAngleParType, ap - mol->par->anglePars);
47274742 }
47284743
47294744 /*
@@ -4737,6 +4752,8 @@
47374752 {
47384753 Molecule *mol;
47394754 Int ival;
4755+ TorsionPar *tp;
4756+ UInt t1, t2, t3, t4;
47404757 Data_Get_Struct(self, Molecule, mol);
47414758 ival = NUM2INT(rb_Integer(val));
47424759 if (ival < -mol->ndihedrals || ival >= mol->ndihedrals)
@@ -4744,7 +4761,14 @@
47444761 if (ival < 0)
47454762 ival += mol->ndihedrals;
47464763 s_RebuildMDParameterIfNecessary(self);
4747- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kDihedralParType, mol->arena->dihedral_par_i[ival]);
4764+ t1 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[ival * 4])->type;
4765+ t2 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[ival * 4 + 1])->type;
4766+ t3 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[ival * 4 + 2])->type;
4767+ t4 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[ival * 4 + 3])->type;
4768+ tp = ParameterLookupDihedralPar(mol->par, t1, t2, t3, t4, 0);
4769+ if (tp == NULL)
4770+ return Qnil;
4771+ return ValueFromMoleculeWithParameterTypeAndIndex(mol, kDihedralParType, tp - mol->par->dihedralPars);
47484772 }
47494773
47504774 /*
@@ -4758,6 +4782,8 @@
47584782 {
47594783 Molecule *mol;
47604784 Int ival;
4785+ TorsionPar *tp;
4786+ UInt t1, t2, t3, t4;
47614787 Data_Get_Struct(self, Molecule, mol);
47624788 ival = NUM2INT(rb_Integer(val));
47634789 if (ival < -mol->nimpropers || ival >= mol->nimpropers)
@@ -4765,7 +4791,14 @@
47654791 if (ival < 0)
47664792 ival += mol->nimpropers;
47674793 s_RebuildMDParameterIfNecessary(self);
4768- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kImproperParType, mol->arena->improper_par_i[ival]);
4794+ t1 = ATOM_AT_INDEX(mol->atoms, mol->impropers[ival * 4])->type;
4795+ t2 = ATOM_AT_INDEX(mol->atoms, mol->impropers[ival * 4 + 1])->type;
4796+ t3 = ATOM_AT_INDEX(mol->atoms, mol->impropers[ival * 4 + 2])->type;
4797+ t4 = ATOM_AT_INDEX(mol->atoms, mol->impropers[ival * 4 + 3])->type;
4798+ tp = ParameterLookupImproperPar(mol->par, t1, t2, t3, t4, 0);
4799+ if (tp == NULL)
4800+ return Qnil;
4801+ return ValueFromMoleculeWithParameterTypeAndIndex(mol, kImproperParType, tp - mol->par->improperPars);
47694802 }
47704803
47714804 /*
@@ -4779,6 +4812,8 @@
47794812 {
47804813 Molecule *mol;
47814814 Int ival;
4815+ VdwPar *vp;
4816+ UInt t1;
47824817 Data_Get_Struct(self, Molecule, mol);
47834818 ival = NUM2INT(rb_Integer(val));
47844819 if (ival < -mol->natoms || ival >= mol->natoms)
@@ -4786,7 +4821,11 @@
47864821 if (ival < 0)
47874822 ival += mol->natoms;
47884823 s_RebuildMDParameterIfNecessary(self);
4789- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kVdwParType, mol->arena->vdw_par_i[ival]);
4824+ t1 = ATOM_AT_INDEX(mol->atoms, ival)->type;
4825+ vp = ParameterLookupVdwPar(mol->par, t1, 0);
4826+ if (vp == NULL)
4827+ return Qnil;
4828+ return ValueFromMoleculeWithParameterTypeAndIndex(mol, kVdwParType, vp - mol->par->vdwPars);
47904829 }
47914830
47924831 /*
Show on old repository browser