计算化学公社

标题: 确定性横向位移分离红细胞，微柱建模程序，不会生成分子模版 [打印本页]

---

作者
Author: zfynpc    时间: 2024-7-23 16:55
标题: 确定性横向位移分离红细胞，微柱建模程序，不会生成分子模版
#include "in_mddefs.h"
void ReadData(char*);
void PrintLAMdata(char*);
void TemplateMaker(char*);

typedef struct {
    VecR r, rv, ra, ract;
    int inChain;
    int beadType;
} Mol;
Mol* mol, * tem;
int nMol, targetFrame, shape_index;
real Rinner, Router; // Diameter or Length of post in Inner or Outer
VecR region, origin;

int main(int argc, char** argv)
{
    int i;
    char bcfg[256], bLAM[256], bTEM[256], shape[256];

    AllocMem(mol, 200000, Mol);
    AllocMem(tem, 200000, Mol);
    targetFrame = 0;

    printf("***********************************************************\n");
    if (argc < 6) {
        printf("5 parameters are needed!\n");
        printf("1. Input File Name; 2. Output File Name; 3. Shape Index; 4. Inner Size; 5. Outer Size\n");
        printf("2D Post Index: 1.SQ; 2.Circle; 3.I-Shape; 4.L-Shape; 5.Inverted L-shape\n");
        exit(1);
    }
    printf("argc = %d\n", argc);
    for (i = 0; i < argc; i++) printf("argv[ %d ] = %s ;\n", i, argv[i]);
    if (argv[1] == NULL) { printf("Need Input Configuration File:\n"); exit(3); }
    else {
        strcpy(bcfg, argv[1]);
        strcpy(bLAM, argv[1]);
        strcat(bcfg, ".dat");
        strcat(bLAM, ".lammpstrj");
    }
    strcpy(bTEM, argv[2]);
    shape_index = strtol(argv[3], NULL, 10);
    Rinner = strtod(argv[4], NULL);
    Router = strtod(argv[5], NULL);

    if (argv[6] == NULL) origin.x = 0.0;
    else origin.x = strtod(argv[6], NULL);
    if (argv[7] == NULL) origin.y = 0.0;
    else origin.y = strtod(argv[7], NULL);
    if (argv[8] == NULL) origin.z = 0.0;
    else origin.z = strtod(argv[8], NULL);

    printf("***********************************************************\n");
    printf("Post Creating !\n");

    ReadData(bcfg);
    TemplateMaker(bTEM);

    free(mol);
    free(tem);
}

void ReadData(char* fileName)
{
    int n, numLine, numPara, paraFlag, countLine, numCFG, numFlag, i, j, j1;
    char s[200];
    char buff[20][20];
    char cfg[30];
    FILE* fp;
    real ra, rb, rx, ry;

    fp = fopen(fileName, "r");
    if (fp == NULL) printf("Error: cant't open the input CFG file!\n");

    numPara = 9;
    countLine = 0;
    numFlag = 0;
    paraFlag = 0;

    do {
        fgets(s, 200, fp);
        if (countLine < 9) {
            sscanf(s, "%s %s %s %s %s", buff[0], buff[1], buff[2], buff[3], buff[4]);
            if (countLine == 1) numCFG = strtol(buff[0], NULL, 10);
            else if (paraFlag == 0 && countLine == 3) nMol = strtol(buff[0], NULL, 10);
            else if (paraFlag == 0 && countLine == 5) {
                ra = strtod(buff[0], NULL);
                rb = strtod(buff[1], NULL);
                region.x = rb - ra;
            }
            else if (paraFlag == 0 && countLine == 6) {
                ra = strtod(buff[0], NULL);
                rb = strtod(buff[1], NULL);
                region.y = rb - ra;
            }
            else if (paraFlag == 0 && countLine == 7) {
                ra = strtod(buff[0], NULL);
                rb = strtod(buff[1], NULL);
                region.z = rb - ra;
            }
            else if (countLine == 8) {
                paraFlag = 1;
                printf("numCFG = %8d;", numCFG);
                printf(" Nmol  = %8d;", nMol);
                printf(" Box.x = %8.4f;", region.x);
                printf(" Box.y = %8.4f", region.y);
                printf(" Box.z = %8.4f\n", region.z);
            }
        }
        if (numCFG == targetFrame) {
            numFlag = 1;
            if (countLine >= 9) {
                sscanf(s, "%s %s %s %s %s", buff[0], buff[1], buff[2], buff[3], buff[4]);
                n = strtol(buff[0], NULL, 10) - 1;
                mol[n].beadType = strtol(buff[1], NULL, 10);
                mol[n].r.x = strtod(buff[2], NULL);
                mol[n].r.y = strtod(buff[3], NULL);
                mol[n].r.z = strtod(buff[4], NULL);
                VWrapAll(mol[n].r);
            }
        }
        ++countLine;
        if (countLine == numPara + nMol) {
            countLine = 0;
            if (numCFG == targetFrame) break;
        }
    } while (!feof(fp));

    if (numFlag == 0) {
        printf("Error: Can't find this frame ***%d***;\n", targetFrame);
        exit(1);
    }

    fclose(fp);
}


void TemplateMaker(char* TEMFile)
{
    int n, NumP, Nnone;
    FILE* fp;
    real ra, rb;

    ra = Router / 2.0;
    rb = Rinner / 2.0;
    NumP = 0;
    Nnone = 0;
    fp = fopen(TEMFile, "w");
    if (fp == NULL) { printf("Error: VMD File Can't Build!\n"); exit(1); }

    for (n = 0; n < nMol; n++) {
        VWrapAll(mol[n].r);
        mol[n].r.x -= origin.x;
        mol[n].r.y -= origin.y;
        mol[n].r.z -= origin.z;
    }

    if (shape_index == 1) {
        //.D Square-shape annulus
        for (n = 0; n < nMol; n++) {
            if (fabs(mol[n].r.x) <= ra && fabs(mol[n].r.y) <= ra) {
                if (fabs(mol[n].r.x) <= rb && fabs(mol[n].r.y) <= rb) {
                    ++Nnone;
                }
                else {
                    VCopy(tem[NumP].r, mol[n].r);
                    tem[NumP].beadType = mol[n].beadType;
                    ++NumP;
                }
            }
        }
    }else if (shape_index == 150) {
        //2D isosceles triangle
        for (n = 0; n < nMol; n++) {
            if (fabs(mol[n].r.x) <= ra && fabs(mol[n].r.y) <= ra) {
                if (3.0 * mol[n].r.y + 4.0 * mol[n].r.x > 15 || 3.0 * mol[n].r.y - 4.0 * mol[n].r.x > 15) {
                    ++Nnone;
                }
                else {
                    VCopy(tem[NumP].r, mol[n].r);
                    tem[NumP].beadType = mol[n].beadType;
                    ++NumP;
                }
            }

        }
    }



    //printf("Shape = %5d\nNumBeads = %5d\nRinner = %.4f\nRouter = %.4f\n", shape_index, NumP, Rinner, Router);
    printf("%d Beads in Post !\n", NumP);

    fprintf(fp, "ITEM: TIMESTEP\n");
    fprintf(fp, "0\n");
    fprintf(fp, "ITEM: NUMBER OF ATOMS\n");
    fprintf(fp, "%d\n", NumP);
    fprintf(fp, "ITEM: BOX BOUNDS\n");
    fprintf(fp, "%8.6f %8.6f\n", -ra, ra);
    fprintf(fp, "%8.6f %8.6f\n", -ra, ra);
    fprintf(fp, "%8.6f %8.6f\n", -0.5 * region.z, 0.5 * region.z);
    fprintf(fp, "ITEM: ATOMS\n");
    for (n = 0; n < NumP; n++) fprintf(fp, "%6d %6d %-8.5f %-8.5f %-8.5f\n", (n + 1), tem[n].beadType, tem[n].r.x, tem[n].r.y, tem[n].r.z);

    fclose(fp);
}

一个用于DLD建模的程序，建立微柱运行程序，提示输入文件读不到，这里的输入文件除了粒子数，盒子大小，还需要什么信息

---

欢迎光临 计算化学公社 (http://bbs.keinsci.com/)

Powered by Discuz! X3.3