Merge pull request #1 from drelatgithub/develop

First runnable program

GraduateProject_DNASelfAssembling/GraduateProject_DNASelfAssembling.cpp

@@ -3,6 +3,8 @@
 
 #include "stdafx.h"
 
+mt19937 gen(123456);
+
 int programInitiation(){
 	ornt2bpornt_init();
 	unitcell_init();
@@ -12,10 +14,8 @@ int _tmain(int argc, _TCHAR* argv[])
 {
 	programInitiation();
 	patchPreset();
-
-
-
+	simulationPrepare();
+	simulationProcess();
 
 	return 0;
 }
-

GraduateProject_DNASelfAssembling/GraduateProject_DNASelfAssembling.vcxproj

@@ -79,15 +79,19 @@
     <Text Include="ReadMe.txt" />
   </ItemGroup>
   <ItemGroup>
+    <ClInclude Include="GraduateProject_Data.h" />
     <ClInclude Include="GraduateProject_PatchesAndPresets.h" />
     <ClInclude Include="GraduateProject_PositionAndOrientation.h" />
+    <ClInclude Include="GraduateProject_Simulation.h" />
     <ClInclude Include="stdafx.h" />
     <ClInclude Include="targetver.h" />
   </ItemGroup>
   <ItemGroup>
+    <ClCompile Include="GraduateProject_Data.cpp" />
     <ClCompile Include="GraduateProject_DNASelfAssembling.cpp" />
     <ClCompile Include="GraduateProject_PatchesAndPresets.cpp" />
     <ClCompile Include="GraduateProject_PositionAndOrientation.cpp" />
+    <ClCompile Include="GraduateProject_Simulation.cpp" />
     <ClCompile Include="stdafx.cpp">
       <PrecompiledHeader Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">Create</PrecompiledHeader>
       <PrecompiledHeader Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">Create</PrecompiledHeader>

GraduateProject_DNASelfAssembling/GraduateProject_DNASelfAssembling.vcxproj.filters

@@ -30,6 +30,12 @@
     <ClInclude Include="GraduateProject_PatchesAndPresets.h">
       <Filter>头文件</Filter>
     </ClInclude>
+    <ClInclude Include="GraduateProject_Simulation.h">
+      <Filter>头文件</Filter>
+    </ClInclude>
+    <ClInclude Include="GraduateProject_Data.h">
+      <Filter>头文件</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="stdafx.cpp">
@@ -44,5 +50,11 @@
     <ClCompile Include="GraduateProject_PatchesAndPresets.cpp">
       <Filter>源文件</Filter>
     </ClCompile>
+    <ClCompile Include="GraduateProject_Simulation.cpp">
+      <Filter>源文件</Filter>
+    </ClCompile>
+    <ClCompile Include="GraduateProject_Data.cpp">
+      <Filter>源文件</Filter>
+    </ClCompile>
   </ItemGroup>
 </Project>

GraduateProject_DNASelfAssembling/GraduateProject_Data.cpp

@@ -0,0 +1,135 @@
+#include"stdafx.h"
+
+const double NN_dH_kcal[4][4][4] = {
+	{ // AX/TY
+		{ 1.2, 2.3, -0.6, -7.6 },
+		{ 5.3, 0.0, -8.4, 0.7 },
+		{ -0.7, -7.8, -3.1, 1.0 },
+		{ -7.2, -1.2, -2.5, -2.7 }
+	},
+	{ // CX/GY
+		{ -0.9, 1.9, -0.7, -8.5 },
+		{ 0.6, -1.5, -8.0, -0.8 },
+		{ -4.0, -10.6, -4.9, -4.1 },
+		{ -7.8, -1.5, -2.8, -5.0 }
+	},
+	{ // GX/CY
+		{ -2.9, 5.2, -0.6, -8.2 },
+		{ -0.7, 3.6, -9.8, 2.3 },
+		{ 0.5, -8.0, -6.0, 3.3 },
+		{ -8.4, 5.2, -4.4, -2.2 }
+	},
+	{ // TX/AY
+		{ 4.7, 3.4, 0.7, -7.2 },
+		{ 7.6, 6.1, -8.2, 1.2 },
+		{ 3.0, -8.5, 1.6, -0.1 },
+		{ -7.6, 1.0, -1.3, 0.2 }
+	}
+};
+const double NN_dS_cal[4][4][4] = {
+	{ // AX/TY
+		{ 1.7, 4.6, -2.3, -21.3 },
+		{ 14.6, -4.4, -22.4, 0.2 },
+		{ -2.3, -21.0, -9.5, 0.9 },
+		{ -20.4, -6.2, -8.3, -10.8 }
+	},
+	{ // CX/GY
+		{ -4.2, 3.7, -2.3, -22.7 },
+		{ -0.6, -7.2, -19.9, -4.5 },
+		{ -13.2, -27.2, -15.3, -11.7 },
+		{ -21.0, -6.1, -8.0, -15.8 }
+	},
+	{ // GX/CY
+		{ -9.8, 14.2, -1.0, -22.2 },
+		{ -3.8, 8.9, -24.4, 5.4 },
+		{ 3.2, -19.9, -15.8, 10.4 },
+		{ -22.4, 13.5, -12.3, -8.4 }
+	},
+	{ // TX/AY
+		{ 12.9, 8.0, 0.7, -21.3 },
+		{ 20.2, 16.4, -22.2, 0.7 },
+		{ 7.4, -22.7, 3.6, -1.7 },
+		{ -21.3, 0.7, -5.3, -1.5 }
+	}
+};
+int NN_dG_check(){
+	double check[4][4][4];
+	int i, j, k;
+	double T37 = 273.15 + 37;
+	for (i = 0; i < 4; i++){
+		for (j = 0; j < 4; j++){
+			for (k = 0; k < 4; k++){
+				check[i][j][k] = NN_dH_kcal[i][j][k] - T37*NN_dS_cal[i][j][k] / 1000;
+				cout << "  \t" << check[i][j][k];
+			}
+			cout << endl;
+		}
+		cout << endl;
+	}
+	return 0;
+}
+const double NN_dH_init_kcal = 0.2;
+const double NN_dS_init_cal = -5.7;
+const double NN_dH_ATPenalty_kcal = 2.2;
+const double NN_dS_ATPenalty_cal = 6.9;
+/*
+	References:
+
+	1. J. SantaLucia, Jr. and D. Hicks, The Thermodynamics of DNA Structural Motifs, Annu. Rev. Biophys. Biomol. Struct. 33, 415 (2004)
+	2. Allawi HT, SantaLucia J Jr. Thermodynamics and NMR of internal G-T mismatches in DNA. Biochemistry 36:10581-94 (1997)
+	3. Allawi HT, SantaLucia J Jr. Nearestneighbor thermodynamics of internal A-C mismatches in DNA: sequence dependence and pH effects. Biochemistry 37:9435-44 (1998)
+	4. Allawi HT, SantaLucia J Jr. Nearestneighbor thermodynamics parameters for internal G-A mismatches in DNA. Biochemistry 37:2170-79 (1998)
+	5. Allawi HT, SantaLucia J Jr. Thermodynamics of internal C-T mismatches in DNA. Nucleic Acids Res. 26:2694-701 (1998)
+	6. Peyret N, Seneviratne PA, Allawi HT, SantaLucia J Jr. Nearest-neighbor thermodynamics and NMR of DNA sequences with internal A-A, C-C, G-G, and T-T mismatches. Biochemistry 38:3468-77 (1999)
+*/
+
+extern const double NN_dH_DE_kcal[2][4][4] = {
+	{
+		{ 0.2, 0.6, -1.1, -6.9 },
+		{ -6.3, -4.4, -5.1, -4.0 },
+		{ -3.7, -4.0, -3.9, -4.9 },
+		{ -2.9, -4.1, -4.2, -0.2 }
+	},
+	{
+		{ -0.5, 4.7, -4.1, -3.8 },
+		{ -5.9, -2.6, -3.2, -5.2 },
+		{ -2.1, -0.2, -3.9, -4.4 },
+		{ -0.7, 4.4, -1.6, 2.9 }
+	}
+};
+extern const double NN_dS_DE_cal[2][4][4] = {
+	{
+		{ 2.3, 3.3, -1.6, -20.0 },
+		{ -17.1, -12.6, -14.0, -10.9 },
+		{ -10.0, -11.9, -10.9, -13.8 },
+		{ -7.6, -13.0, -15.0, -0.5 }
+	},
+	{
+		{ -1.1, 14.2, -13.1, -12.6 },
+		{ -16.5, -7.4, -10.4, -15.0 },
+		{ -3.9, -0.1, -11.2, -13.1 },
+		{ -0.8, 14.9, -3.6, 10.4 }
+	}
+};
+int NN_dG_DE_check(){
+	double check[2][4][4];
+	int i, j, k;
+	double T37 = 273.15 + 37;
+	for (i = 0; i < 2; i++){
+		for (j = 0; j < 4; j++){
+			for (k = 0; k < 4; k++){
+				check[i][j][k] = NN_dH_DE_kcal[i][j][k] - T37*NN_dS_DE_cal[i][j][k] / 1000;
+				cout << "  \t" << check[i][j][k];
+			}
+			cout << endl;
+		}
+		cout << endl;
+	}
+	return 0;
+}
+/*
+	References:
+
+	1. J. SantaLucia, Jr. and D. Hicks, The Thermodynamics of DNA Structural Motifs, Annu. Rev. Biophys. Biomol. Struct. 33, 415 (2004)
+	2. Bommarito S, Peyret N, SantaLucia J Jr. Thermodynamic parameters for DNA sequences with dangling ends. Nucleic Acids Res. 28:1929-34 (2000)
+*/

GraduateProject_DNASelfAssembling/GraduateProject_Data.h

@@ -0,0 +1,38 @@
+#pragma once
+
+extern const double NN_dH_kcal[4][4][4];
+extern const double NN_dS_cal[4][4][4];
+/*
+	NN_dH_kcal[a][b][c], NN_dS_cal[a][b][c]
+
+	5'  ---  a  ---  b  ---  3'
+	         |       |
+	3'  ---  a' ---  c  ---  5'
+
+	where a' is the paired nucleotide of a.
+	If c = b' then the index [a][b][b'] and [b'][a'][a] must refer to the same value.
+*/
+//int NN_dG_check();
+extern const double NN_dH_init_kcal;
+extern const double NN_dS_init_cal;
+extern const double NN_dH_ATPenalty_kcal;
+extern const double NN_dS_ATPenalty_cal;
+
+extern const double NN_dH_DE_kcal[2][4][4];
+extern const double NN_dS_DE_cal[2][4][4];
+/*
+	NN_dH_DE_kcal[0,1][a][b], NN_dS_DE_cal[0,1][a][b]
+
+	0 mode: 5' dangling ends
+
+		5'  ---  b  ---  a  ---  3'
+		                 |
+						 a' ---  5'
+
+	1 mode: 3' dangling ends
+
+		                 a' ---  3'
+						 |
+		3'  ---  b  ---  a  ---  5'
+*/
+//int NN_dG_DE_check();

GraduateProject_DNASelfAssembling/GraduateProject_PatchesAndPresets.cpp

@@ -5,10 +5,7 @@ atom_in_unitcell atom_in_unitcell::set(int nx, int ny, int nz, int nsd, int norn
 	sd = nsd;
 	int mainOrnt = nornt0;
 	int subOrnt = 0, subOrnt_as_pow = 7 - nornt0^nornt1;
-	while (subOrnt_as_pow > 1){
-		subOrnt_as_pow >>= 1;
-		subOrnt++;
-	}
+	for (; subOrnt_as_pow > 1; subOrnt_as_pow >>= 1)subOrnt++;
 	ornt = mainOrnt + 8 * subOrnt;
 	return *this;
 }
@@ -28,8 +25,10 @@ int unitcell_init(){
 DNAmol::DNAmol(){
 	px.x = px.y = px.z = 0;
 	ornt = 0;
-	for (int i = 0; i < 4; i++){
+	int i, j;
+	for (i = 0; i < 4; i++){
 		correctbond[i] = -1;
+		for (j = 0; j < 8; j++)patch[i][j] = 'T';
 	}
 }
 DNAmol DNAmol::put(int nx, int ny, int nz, int nornt){
@@ -42,10 +41,53 @@ DNAmol DNAmol::put(const atom_in_unitcell &a, int x_init, int y_init, int z_init
 	ornt = a.ornt;
 	return *this;
 }
+int DNAmol::displayPatch()const{
+	for (int i = 0; i < 4; i++){
+		cout.write(patch[i], 8);
+	}
+	cout << endl;
+	return 0;
+}
+int DNAmol::displayPatch(ofstream &out)const{
+	for (int i = 0; i < 4; i++){
+		out.write(patch[i], 8);
+	}
+	out << endl;
+	return 0;
+}
+int DNAmol::findPatchSerial(int whichOrnt)const{ // Find the patch serial in the designated orientation. -1 if not found.
+	for (int m = 0; m < 4; m++){
+		if (ornt2bpornt[ornt][m] == whichOrnt){
+			return m;
+		}
+	}
+	return -1;
+}
 int stage[_Nx][_Ny][_Nz]; // To store DNAmol serial. -1 if not occupied. Using periodic boundary conditions.
 DNAmol *mol;
 int N; // Total molecules.
 
+char nt[4] = { 'A', 'C', 'G', 'T' };
+int ntSerial(char which_nt){
+	switch (which_nt){
+	case'A':return 0;
+	case'C':return 1;
+	case'G':return 2;
+	case'T':return 3;
+	}
+}
+int ntSerialPair[4] = { 3, 2, 1, 0 };
+char ntPair(char which_nt){
+	return nt[ntSerialPair[ntSerial(which_nt)]];
+}
+int couldPatchInteract[4][4] = {
+	{ 0, 1, 1, 0 },
+	{ 1, 0, 0, 1 },
+	{ 1, 0, 0, 1 },
+	{ 0, 1, 1, 0 }
+};
+int isPatchConsequent[4] = { 1, 0, 0, 1 };
+
 int stageClear(){
 	int i, j, k;
 	for (i = 0; i < _Nx; i++){
@@ -89,22 +131,19 @@ int correctbonding(int n0serial, int n01x, int n01y, int n01z){
 	int ornt0 = 4 * (n01x + 1) / 2 + 2 * (n01y + 1) / 2 + (n01z + 1) / 2;
 	int ornt1 = 7 - ornt0;
 	int n0ps = -1, n1ps = -1; // patch serial for the right bonding. -1 if not exist.
-	int m;
-	for (m = 0; m < 4; m++){
-		if (ornt2bpornt[mol[n0serial].ornt][m] == ornt0){
-			n0ps = m; break;
-		}
-	}
-	for (m = 0; m < 4; m++){
-		if (ornt2bpornt[mol[n1serial].ornt][m] == ornt1){
-			n1ps = m; break;
-		}
-	}
-	if (n0ps >= 0 && n1ps >= 0 && abs(n0ps - n1ps) == 2){
+	int temp_ntSerial;
+	n0ps = mol[n0serial].findPatchSerial(ornt0);
+	n1ps = mol[n1serial].findPatchSerial(ornt1);
+	if (n0ps >= 0 && n1ps >= 0 && couldPatchInteract[n0ps][n1ps]){
 		mol[n0serial].correctbond[n0ps] = n1serial;
 		mol[n1serial].correctbond[n1ps] = n0serial;
 		// patch assign
-
+		static uniform_int_distribution<> patchDis(0, 3);
+		for (int i = 0; i < 8; i++){
+			temp_ntSerial = patchDis(gen);
+			mol[n0serial].patch[n0ps][i] = nt[temp_ntSerial];
+			mol[n1serial].patch[n1ps][7 - i] = nt[3 - temp_ntSerial];
+		}
 		return 0;
 	}
 	return 1;
@@ -124,11 +163,17 @@ int nucleotidePreset(){
 			}
 		}
 	}
+	ofstream out("F:\\DNAPatches.txt");
+	for (m = 0; m < N; m++){
+		mol[m].displayPatch(out);
+	}
+	out.close();
 	return 0;
 }
 int patchPreset(){
 	stageClear();
 	molPreset();
 	nucleotidePreset();
+	stageClear();
 	return 0;
 }

GraduateProject_DNASelfAssembling/GraduateProject_PatchesAndPresets.h

@@ -17,9 +17,19 @@ struct DNAmol{
 	DNAmol();
 	DNAmol put(int nx, int ny, int nz, int nornt);
 	DNAmol put(const atom_in_unitcell &a, int x_init, int y_init, int z_init);
+	int displayPatch()const;
+	int displayPatch(ofstream &out)const;
+	int findPatchSerial(int whichOrnt)const; // Find the patch serial in the designated orientation. -1 if not found.
 };
 extern DNAmol *mol;
 extern int stage[_Nx][_Ny][_Nz]; // To store DNAmol serial. -1 if not occupied. Using periodic boundary conditions.
 extern int N; // Total molecules.
 
+extern char nt[4];
+int ntSerial(char which_nt);
+extern int ntSerialPair[4];
+char ntPair(char which_nt);
+extern int couldPatchInteract[4][4];
+extern int isPatchConsequent[4]; // If the nucleotide sequence in a patch is 5'-3', it is consequent.
+
 int patchPreset();