Search in sources :

Example 1 with SharedDoubleArray

use of edu.rit.pj.reduction.SharedDoubleArray in project ffx by mjschnie.

the class GeneralizedKirkwood method initAtomArrays.

private void initAtomArrays() {
    if (fixedRadii) {
        fixedRadii = false;
    }
    x = particleMeshEwald.coordinates[0][0];
    y = particleMeshEwald.coordinates[0][1];
    z = particleMeshEwald.coordinates[0][2];
    globalMultipole = particleMeshEwald.globalMultipole[0];
    inducedDipole = particleMeshEwald.inducedDipole[0];
    inducedDipoleCR = particleMeshEwald.inducedDipoleCR[0];
    neighborLists = particleMeshEwald.neighborLists;
    grad = particleMeshEwald.getGradient();
    torque = particleMeshEwald.getTorque();
    lambdaGrad = particleMeshEwald.getLambdaGradient();
    lambdaTorque = particleMeshEwald.getLambdaTorque();
    if (sharedGKField[0] == null || sharedGKField[0].length() < nAtoms) {
        sharedGKField[0] = new SharedDoubleArray(nAtoms);
        sharedGKField[1] = new SharedDoubleArray(nAtoms);
        sharedGKField[2] = new SharedDoubleArray(nAtoms);
        sharedGKFieldCR[0] = new SharedDoubleArray(nAtoms);
        sharedGKFieldCR[1] = new SharedDoubleArray(nAtoms);
        sharedGKFieldCR[2] = new SharedDoubleArray(nAtoms);
        sharedBornGrad = new SharedDoubleArray(nAtoms);
        baseRadius = new double[nAtoms];
        baseRadiusWithBondi = new double[nAtoms];
        overlapScale = new double[nAtoms];
        rDisp = new double[nAtoms];
        born = new double[nAtoms];
        use = new boolean[nAtoms];
    }
    fill(use, true);
    for (int i = 0; i < nAtoms; i++) {
        baseRadius[i] = 2.0;
        // Original value based on small molecule parameterization.
        overlapScale[i] = DEFAULT_OVERLAP_SCALE;
        // overlapScale[i] = 0.60;     // New default value based on 2016 amino acid GK parameterization.
        if (useFittedRadii) {
            overlapScale[i] = solventRadii.getOverlapScale();
        }
        int atomicNumber = atoms[i].getAtomicNumber();
        AtomType atomType = atoms[i].getAtomType();
        switch(atomicNumber) {
            case 0:
                baseRadius[i] = 0.0;
                break;
            case 1:
                baseRadius[i] = 1.2;
                overlapScale[i] = hydrogenOverlapScale;
                break;
            case 2:
                baseRadius[i] = 1.4;
                break;
            case 5:
                baseRadius[i] = 1.8;
                break;
            case 6:
                baseRadius[i] = 1.7;
                break;
            case 7:
                baseRadius[i] = 1.55;
                break;
            case 8:
                baseRadius[i] = 1.52;
                break;
            case 9:
                baseRadius[i] = 1.47;
                break;
            case 10:
                baseRadius[i] = 1.54;
                break;
            case 14:
                baseRadius[i] = 2.1;
                break;
            case 15:
                baseRadius[i] = 1.8;
                break;
            case 16:
                baseRadius[i] = 1.8;
                break;
            case 17:
                baseRadius[i] = 1.75;
                break;
            case 18:
                baseRadius[i] = 1.88;
                break;
            case 34:
                baseRadius[i] = 1.9;
                break;
            case 35:
                baseRadius[i] = 1.85;
                break;
            case 36:
                baseRadius[i] = 2.02;
                break;
            case 53:
                baseRadius[i] = 1.98;
                break;
            case 54:
                baseRadius[i] = 2.16;
                break;
            default:
                baseRadius[i] = 2.00;
        }
        double bondiFactor = bondiScale;
        int atomNumber = atoms[i].getIndex() + 1;
        if (useFittedRadii) {
            // First check to see if this atom is in the hardcoded maps.
            switch(radiiMapType) {
                default:
                case ATOMTYPE:
                    // Check for hard-coded AtomType bondi factor.
                    if (solventRadii.getAtomBondiMap().containsKey(atomType.type)) {
                        bondiFactor = solventRadii.getAtomBondiMap().get(atomType.type);
                        if (verboseRadii) {
                            logger.info(String.format(" (GK) TypeToBondi: Atom %3s-%-4s (%d) with AtomType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
                        }
                    }
                    // Many aminos use types 8,12 for their CA,HA so these can't be specified in the map.
                    if ((atoms[i].getAtomType().type == 8 || atoms[i].getAtomType().type == 12)) {
                        if (atoms[i].getResidueName().equals("ALA")) {
                            bondiFactor = 1.60;
                            if (verboseRadii) {
                                logger.info(String.format(" (GK) TypeToBondi: Atom %3s-%-4s (%d) with AtomType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
                            }
                        }
                    }
                    // Currently, we only want S+HS on CYS.  So the CB,HB (types 43,44) are disabled in the table.
                    if ((atoms[i].getAtomType().type == 43 || atoms[i].getAtomType().type == 44)) {
                        if (atoms[i].getResidueName().equals("CYD")) {
                            bondiFactor = 1.02;
                            if (verboseRadii) {
                                logger.info(String.format(" (GK) TypeToBondi: Atom %3s-%-4s (%d) with AtomType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
                            }
                        }
                    }
                    break;
                case BIOTYPE:
                    Map<String, BioType> bioTypes = forceField.getBioTypeMap();
                    BioType bioType = null;
                    for (BioType one : bioTypes.values()) {
                        if (one.atomType == atomType.type) {
                            bioType = one;
                            break;
                        }
                    }
                    if (bioType == null) {
                        logger.severe(String.format("Couldn't find biotype for %s", atomType, toString()));
                    }
                    // Check for hard-coded BioType bondi factor.
                    if (solventRadii.getBioBondiMap().containsKey(bioType.index)) {
                        double factor = solventRadii.getBioBondiMap().get(bioType.index);
                        bondiFactor = factor;
                        if (verboseRadii) {
                            logger.info(String.format(" (GK) BiotypeToBondi: Atom %3s-%-4s (%d) with BioType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, bioType.index, baseRadius[i] * bondiFactor, bondiFactor));
                        }
                    }
                    // Many aminos use types 8,12 for their CA,HA so these can't be specified in the map.
                    if (bioType.index == 8 || bioType.index == 12) {
                        if (atoms[i].getResidueName().equals("ALA")) {
                            bondiFactor = 1.60;
                            if (verboseRadii) {
                                logger.info(String.format(" (GK) BiotypeToBondi: Atom %3s-%-4s (%d) with BioType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, bioType.index, baseRadius[i] * bondiFactor, bondiFactor));
                            }
                        }
                    }
                    // Currently, we only want S+HS on CYS.  So the CB,HB (types 43,44) are disabled in the table.
                    if (bioType.index == 83 || bioType.index == 84) {
                        if (atoms[i].getResidueName().equals("CYD")) {
                            bondiFactor = 1.02;
                            if (verboseRadii) {
                                logger.info(String.format(" (GK) BiotypeToBondi: Atom %3s-%-4s (%d) with BioType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, bioType.index, baseRadius[i] * bondiFactor, bondiFactor));
                            }
                        }
                    }
                    break;
            }
        }
        // Next, check if this Atom has an ISolvRad forcefield parameter.
        // if (atoms[i].getISolvRadType() != null) {
        // bondiFactor = atoms[i].getISolvRadType().radiusScale;
        // logger.info(String.format(" (GK) ISolvRad parameter for Atom %3s-%-4s with AtomType %d to %.2f (bondi factor %.2f)",
        // atoms[i].getResidueName(), atoms[i].getName(), atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
        // }
        ISolvRadType iSolvRadType = forceField.getISolvRadType(Integer.toString(atomType.type));
        if (iSolvRadType != null) {
            bondiFactor = iSolvRadType.radiusScale;
            if (verboseRadii) {
                logger.info(String.format(" (GK) ISolvRad parameter for Atom %3s-%-4s with AtomType %d to %.2f (bondi factor %.2f)", atoms[i].getResidueName(), atoms[i].getName(), atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
            }
        }
        // Finally, check for command-line bondi factor override.
        if (radiiOverride.containsKey(atomType.type) && !radiiByNumberMap.containsKey(atomNumber)) {
            bondiFactor = radiiOverride.get(atomType.type);
            logger.info(String.format(" (GK) Scaling Atom %3s-%-4s with AtomType %d to %.2f (bondi factor %.2f)", atoms[i].getResidueName(), atoms[i].getName(), atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
        }
        if (radiiByNumberMap.containsKey(atomNumber)) {
            bondiFactor = radiiByNumberMap.get(atomNumber);
            logger.info(String.format(" (GK) Scaling Atom number %d, %3s-%-4s, with factor %.2f", atomNumber, atoms[i].getResidueName(), atoms[i].getName(), bondiFactor));
        }
        if (!radiiOverride.containsKey(atomType.type)) {
            bondiFactor *= globalRadiiScale;
        } else {
            logger.fine(String.format(" Not scaling atom type %d", atomType.type));
        }
        baseRadiusWithBondi[i] = baseRadius[i] * bondiFactor;
    }
    if (dispersionRegion != null) {
        dispersionRegion.init();
    }
    if (cavitationRegion != null) {
        cavitationRegion.init();
    }
}
Also used : SharedDoubleArray(edu.rit.pj.reduction.SharedDoubleArray) AtomType(ffx.potential.parameters.AtomType) ISolvRadType(ffx.potential.parameters.ISolvRadType) BioType(ffx.potential.parameters.BioType)

Aggregations

SharedDoubleArray (edu.rit.pj.reduction.SharedDoubleArray)1 AtomType (ffx.potential.parameters.AtomType)1 BioType (ffx.potential.parameters.BioType)1 ISolvRadType (ffx.potential.parameters.ISolvRadType)1