__init__.py 59 KB
Newer Older
1001
            os.makedirs(d)
1002
        self._restraint_filename = "%s/%s.restraint.txt" % (d, prefix)
1003
1004
1005
1006
        self._bond_filename = "%s/%s.bonds.txt" % (d, prefix)
        self._angle_filename = "%s/%s.angles.txt" % (d, prefix)
        self._dihedral_filename = "%s/%s.dihedrals.txt" % (d, prefix)
        self._exclusion_filename = "%s/%s.exculsions.txt" % (d, prefix)
1007
        self._bond_angle_filename = "%s/%s.bondangles.txt" % (d, prefix)
1008
        self._product_potential_filename = "%s/%s.prodpot.txt" % (d, prefix)
cmaffeo2's avatar
cmaffeo2 committed
1009
        self._group_sites_filename = "%s/%s.groups.txt" % (d, prefix)
1010
        
1011
1012
        # self._writeArbdCoordFile( prefix + ".coord.txt" )
        self._writeArbdParticleFile( prefix + ".particles.txt" )
1013
        self._writeArbdRestraintFile()
1014
1015
1016
1017
        self._writeArbdBondFile()
        self._writeArbdAngleFile()
        self._writeArbdDihedralFile()
        self._writeArbdExclusionFile()
1018
        self._writeArbdBondAngleFile()
1019
        self._writeArbdProductPotentialFile()
cmaffeo2's avatar
cmaffeo2 committed
1020
        self._writeArbdGroupSitesFile()
1021
        self._writeArbdPotentialFiles( prefix, directory = d )
1022
        self._writeArbdConf( prefix, numSteps=numSteps, outputPeriod=outputPeriod, restart_file=restart_file )
1023
        
1024
1025
1026
1027
1028
1029
    # def _writeArbdCoordFile(self, filename):
    #     with open(filename,'w') as fh:
    #         for p in self.particles:
    #             fh.write("%f %f %f\n" % tuple(x for x in p.collapsedPosition()))

    def _writeArbdParticleFile(self, filename):
1030
        with open(filename,'w') as fh:
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
            if self.particle_integrator == "Brown":
                for p in self.particles:
                    data = tuple([p.idx,p.type_.name] + [x for x in p.collapsedPosition()])
                    fh.write("ATOM %d %s %f %f %f\n" % data)
            else:
                for p in self.particles:
                    data = [p.idx,p.type_.name] + [x for x in p.collapsedPosition()]
                    try:
                        data = data + p.momentum
                    except:
                        try:
                            data = data + p.velocity*p.mass
                        except:
                            data = data + [0,0,0]
                    fh.write("ATOM %d %s %f %f %f %f %f %f\n" % tuple(data))
                
    def _write_rigid_group_file(self, filename, groups):
        with open(filename,'w') as fh:
            for g in groups:
                fh.write("#Group\n")
                try:
                    if len(g.trans_damping) != 3: raise
                    fh.write(" ".join(str(v) for v in g.trans_damping) + " ")
                except:
                    raise ValueError("Group {} lacks 3-value 'trans_damping' attribute")
                try:
                    if len(g.rot_damping) != 3: raise
                    fh.write(" ".join(str(v) for v in g.rot_damping) + " ")
                except:
                    raise ValueError("Group {} lacks 3-value 'rot_damping' attribute")
                fh.write("{}\n".format(len(g)))
                particles = [p for p in g]

                def chunks(l, n):
                    """Yield successive n-sized chunks from l."""
                    for i in range(0, len(l), n):
                        yield l[i:i + n]

                for c in chunks(particles,8):
                    fh.write(" ".join(str(p.idx) for p in c) + "\n")


    def _writeArbdConf(self, prefix, randomSeed=None, numSteps=100000000, outputPeriod=10000, restart_file=None):
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
        ## TODO: raise exception if _writeArbdPotentialFiles has not been called
        filename = "%s.bd" % prefix

        ## Prepare a dictionary to fill in placeholders in the configuration file
        params = self.__dict__.copy() # get parameters from System object

        if randomSeed is None:
            params['randomSeed']     = ""
        else:
            params['randomSeed'] = "seed %s" % randomSeed
        params['numSteps']       = int(numSteps)

1086
1087
        # params['coordinateFile'] = "%s.coord.txt" % prefix
        params['particleFile'] = "%s.particles.txt" % prefix
1088
        if restart_file is None:
1089
1090
            params['restartCoordinates'] = ""
        else:
1091
            params['restartCoordinates'] = "restartCoordinates %s" % restart_file
1092
1093
1094
1095
        params['outputPeriod'] = outputPeriod

        for k,v in zip('XYZ', self.dimensions):
            params['dim'+k] = v
1096
1097
1098
1099
1100
1101
1102

        if self.origin is None:
            for k,v in zip('XYZ', self.dimensions):
                params['origin'+k] = -v*0.5
        else:
            for k,v in zip('XYZ', self.origin):
                params['origin'+k] = v
cmaffeo2's avatar
cmaffeo2 committed
1103
        
1104
        params['pairlist_distance'] -= params['cutoff'] 
1105

1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
        """ Find rigid groups """
        rigid_groups = []
        def get_rigid_groups(parent):
            ret_list = []
            for c in parent.children:
                is_rigid = c.rigid if 'rigid' in c.__dict__ else False
                if is_rigid:
                    rigid_groups.append(c)
                elif isinstance(c,Group):
                    get_rigid_groups(c)
        get_rigid_groups(self)

        if len(rigid_groups) > 0:
            self.particle_integrator = 'FusDynamic'
            rb_group_filename = "{}.rb-group.txt".format(prefix)
            params['particle_integrator'] = """FusDynamics
groupFileName {}
scaleFactor 0.05""".format(rb_group_filename)
            self._write_rigid_group_file(rb_group_filename, rigid_groups)

1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
        ## Actually write the file
        with open(filename,'w') as fh:
            fh.write("""{randomSeed}
timestep {timestep}
steps {numSteps}
numberFluct 0                   # deprecated

interparticleForce 1            # other values deprecated
fullLongRange 0                 # deprecated
temperature {temperature}
1136
ParticleDynamicType {particle_integrator}
1137
1138
1139
1140
1141
1142
1143

outputPeriod {outputPeriod}
## Energy doesn't actually get printed!
outputEnergyPeriod {outputPeriod}
outputFormat dcd

## Infrequent domain decomposition because this kernel is still very slow
1144
decompPeriod {decomp_period}
1145
cutoff {cutoff}
1146
pairlistDistance {pairlist_distance}
1147
1148
1149

origin {originX} {originY} {originZ}
systemSize {dimX} {dimY} {dimZ}
1150
1151

{extra_bd_file_lines}
1152
1153
1154
1155
\n""".format(**params))
            
            ## Write entries for each type of particle
            for pt,num in self.getParticleTypesAndCounts():
cmaffeo2's avatar
cmaffeo2 committed
1156
                ## TODO create new particle types if existing has grid
1157
1158
                particleParams = pt.__dict__.copy()
                particleParams['num'] = num
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
                if self.particle_integrator in ('Brown','Brownian'):
                    try:
                        D = pt.diffusivity
                    except:
                        """ units "k K/(amu/ns)" "AA**2/ns" """
                        D = 831447.2 * self.temperature / (pt.mass * pt.damping_coefficient)
                    particleParams['dynamics'] = 'diffusion {D}'.format(D = D)
                elif self.particle_integrator in ('Langevin','FusDynamic'):
                    try:
                        gamma = pt.damping_coefficient
                    except:
                        """ units "k K/(AA**2/ns)" "amu/ns" """
                        gamma = 831447.2 * self.temperature / (pt.mass*pt.diffusivity)
                    particleParams['dynamics'] = """mass {mass}
transDamping {g} {g} {g}
""".format(mass=pt.mass, g=gamma)
                else:
                    raise ValueError("Unrecognized particle integrator '{}'".format(self.particle_integrator))
1177
1178
1179
                fh.write("""
particle {name}
num {num}
1180
{dynamics}
1181
1182
""".format(**particleParams))
                if 'grid' in particleParams:
cmaffeo2's avatar
cmaffeo2 committed
1183
                    if not isinstance(pt.grid, list): pt.grid = [pt.grid]
1184
1185
1186
1187
1188
1189
1190
1191
                    for g,s in pt.grid:
                        ## TODO, use Path.relative_to?
                        try:
                            fh.write("gridFile {}\n".format(g.relative_to(os.getcwd())))
                        except:
                            fh.write("gridFile {}\n".format(g))

                        fh.write("gridFileScale {}\n".format(s))
cmaffeo2's avatar
cmaffeo2 committed
1192

1193
                else:
1194
                    fh.write("gridFile {}/null.dx\n".format(self.potential_directory))
1195
1196
1197
1198

            ## Write coordinates and interactions
            fh.write("""
## Input coordinates
1199
inputParticles {particleFile}
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
{restartCoordinates}

## Interaction potentials
tabulatedPotential  1
## The i@j@file syntax means particle type i will have NB interactions with particle type j using the potential in file
""".format(**params))
            for pair,f in zip(self._particleTypePairIter(), self._nbParamFiles):
                i,j,t1,t2 = pair
                fh.write("tabulatedFile %d@%d@%s\n" % (i,j,f))

            ## Bonded interactions
1211
            restraints = self.get_restraints()
1212
1213
1214
1215
            bonds = self.get_bonds()
            angles = self.get_angles()
            dihedrals = self.get_dihedrals()
            exclusions = self.get_exclusions()
1216
            bond_angles = self.get_bond_angles()
1217
            prod_pots = self.get_product_potentials()
cmaffeo2's avatar
cmaffeo2 committed
1218
1219
1220
            # group_sites = self.get_group_sites()
            group_sites = self.group_sites

1221
            if len(bonds) > 0:
1222
                for b in self._get_bond_potentials():
1223
1224
1225
                    fh.write("tabulatedBondFile %s\n" % b)

            if len(angles) > 0:
1226
                for b in self._get_angle_potentials():
1227
1228
1229
1230
1231
1232
                    fh.write("tabulatedAngleFile %s\n" % b)

            if len(dihedrals) > 0:
                for b in list(set([b for i,j,k,l,b in dihedrals])):
                    fh.write("tabulatedDihedralFile %s\n" % b)

1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
            if len(restraints) > 0:
                fh.write("inputRestraints %s\n" % self._restraint_filename)
            if len(bonds) > 0:
                fh.write("inputBonds %s\n" % self._bond_filename)
            if len(angles) > 0:
                fh.write("inputAngles %s\n" % self._angle_filename)
            if len(dihedrals) > 0:
                fh.write("inputDihedrals %s\n" % self._dihedral_filename)
            if len(exclusions) > 0:
                fh.write("inputExcludes %s\n" % self._exclusion_filename)
1243
1244
            if len(bond_angles) > 0:
                fh.write("inputBondAngles %s\n" % self._bond_angle_filename)
1245
1246
            if len(prod_pots) > 0:
                fh.write("inputProductPotentials %s\n" % self._product_potential_filename)
cmaffeo2's avatar
cmaffeo2 committed
1247
1248
            if len(group_sites) > 0:
                fh.write("inputGroups %s\n" % self._group_sites_filename)
1249
1250
1251
1252
     
        write_null_dx = False
        for pt,num in self.getParticleTypesAndCounts():
            if "grid" not in pt.__dict__: 
1253
1254
                gridfile = "{}/null.dx".format(self.potential_directory)
                with open(gridfile, 'w') as fh:
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
                    fh.write("""object 1 class gridpositions counts  2 2 2
origin {originX} {originY} {originZ}
delta  {dimX} 0.000000 0.000000
delta  0.000000 {dimY} 0.000000
delta  0.000000 0.000000 {dimZ}
object 2 class gridconnections counts  2 2 2
object 3 class array type float rank 0 items 8 data follows
0.0	0.0	0.0	
0.0	0.0	0.0	
0.0	0.0	
attribute "dep" string "positions"
object "density" class field 
component "positions" value 1
component "connections" value 2
component "data" value 3
""".format(**params))
                    break

    def getParticleTypesAndCounts(self):
        ## TODO: remove(?)
        return sorted( self.type_counts.items(), key=lambda x: x[0] )

    def _particleTypePairIter(self):
        typesAndCounts = self.getParticleTypesAndCounts()
        for i in range(len(typesAndCounts)):
            t1 = typesAndCounts[i][0]
            for j in range(i,len(typesAndCounts)):
                t2 = typesAndCounts[j][0]
                yield( (i,j,t1,t2) )
    
    def _writeArbdPotentialFiles(self, prefix, directory = "potentials"):
        try: 
            os.makedirs(directory)
        except OSError:
            if not os.path.isdir(directory):
                raise

        pathPrefix = "%s/%s" % (directory,prefix)
        self._writeNonbondedParameterFiles( pathPrefix + "-nb" )
        # self._writeBondParameterFiles( pathPrefix )
        # self._writeAngleParameterFiles( pathPrefix )
        # self._writeDihedralParameterFiles( pathPrefix )
                
    def _writeNonbondedParameterFiles(self, prefix):
        x = np.arange(0, self.cutoff, self.nbResolution)
        for i,j,t1,t2 in self._particleTypePairIter():
            f = "%s.%s-%s.dat" % (prefix, t1.name, t2.name)
            scheme = self._getNbScheme(t1,t2)
            scheme.write_file(f, t1, t2, rMax = self.cutoff)
            self._nbParamFiles.append(f)

    def _getNonbondedPotential(self,x,a,b):
        return a*(np.exp(-x/b))    

1309
1310
1311
1312
1313
1314
    def _writeArbdRestraintFile( self ):
        with open(self._restraint_filename,'w') as fh:
            for i,restraint in self.get_restraints():
                item = [i.idx]
                if len(restraint) == 1:
                    item.append(restraint[0])
1315
                    item.extend(i.collapsedPosition())
1316
1317
1318
1319
1320
1321
1322
                elif len(restraint) == 2:
                    item.append(restraint[0])
                    item.extend(restraint[1])
                elif len(restraint) == 5:
                    item.extend(restraint)
                fh.write("RESTRAINT %d %f %f %f %f\n" % tuple(item))

1323
    def _writeArbdBondFile( self ):
1324
        for b in self._get_bond_potentials():
1325
            if type(b) is not str and not isinstance(b, Path):
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
                b.write_file()

        with open(self._bond_filename,'w') as fh:
            for i,j,b,ex in self.get_bonds():
                item = (i.idx, j.idx, str(b))
                if ex:
                    fh.write("BOND REPLACE %d %d %s\n" % item)
                else:
                    fh.write("BOND ADD %d %d %s\n" % item)

    def _writeArbdAngleFile( self ):
1337
        for b in self._get_angle_potentials():
1338
            if type(b) is not str and not isinstance(b, Path):
1339
1340
1341
1342
1343
1344
1345
1346
1347
                b.write_file()

        with open(self._angle_filename,'w') as fh:
            for b in self.get_angles():
                item = tuple([p.idx for p in b[:-1]] + [str(b[-1])])
                fh.write("ANGLE %d %d %d %s\n" % item)

    def _writeArbdDihedralFile( self ):
        for b in list( set( [b for i,j,k,l,b in self.get_dihedrals()] ) ):
1348
            if type(b) is not str and not isinstance(b, Path):
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
                b.write_file()

        with open(self._dihedral_filename,'w') as fh:
            for b in self.get_dihedrals():
                item = tuple([p.idx for p in b[:-1]] + [str(b[-1])])
                fh.write("DIHEDRAL %d %d %d %d %s\n" % item)

    def _writeArbdExclusionFile( self ):
        with open(self._exclusion_filename,'w') as fh:
            for ex in self.get_exclusions():
                item = tuple(int(p.idx) for p in ex)
                fh.write("EXCLUDE %d %d\n" % item)
cmaffeo2's avatar
cmaffeo2 committed
1361

1362
1363
1364
1365
1366
    def _writeArbdBondAngleFile( self ):
        if len(self.bond_angles) > 0:
            with open(self._bond_angle_filename,'w') as fh:
                for b in self.get_bond_angles():
                    item = tuple([p.idx for p in b[:-1]] + [str(p) for p in b[-1]])
1367
                    fh.write("BONDANGLE %d %d %d %d %s %s %s\n" % item)
1368

1369
1370
1371
    def _writeArbdProductPotentialFile( self ):
        if len(self.product_potentials) > 0:
            with open(self._product_potential_filename,'w') as fh:
cmaffeo2's avatar
cmaffeo2 committed
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
                for pot in self.get_product_potentials():
                    line = "PRODUCTPOTENTIAL "
                    for ijk_tb in pot:
                        ijk = ijk_tb[:-1]
                        tb = ijk_tb[-1]
                        if type(tb) is tuple or type(tb) is list:
                            if len(tb) != 2: raise ValueError("Invalid product potential")
                            type_,b = tb
                            if type(type_) is not str: raise ValueError("Invalid product potential: unrecognized specification of potential type")
                        else:
                            type_ = ""
                            b = tb
                        if type(b) is not str and not isinstance(b, Path):
                            b.write_file()
                        line = line+" ".join([str(x.idx) for x in ijk])+" "
                        line = line+" ".join([str(x) for x in [type_,b] if x != ""])+" "
                    fh.write(line)
1389

cmaffeo2's avatar
cmaffeo2 committed
1390
1391
1392
1393
1394
1395
1396
    def _writeArbdGroupSitesFile( self ):
        if len(self.group_sites) > 0:
            with open(self._group_sites_filename,'w') as fh:
                for i,g in enumerate(self.group_sites):
                    assert( i+len(self.particles) == g.idx )
                    ids = " ".join([str(int(p.idx)) for p in g.particles])
                    fh.write("GROUP %s\n" % ids)
1397

1398
    def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
1399
1400
1401
        ## TODO: cache coordinates using numpy arrays for quick min/max
        raise(NotImplementedError)

1402
    def write_namd_configuration( self, output_name, minimization_steps=4800, num_steps = 1e6,
1403
1404
1405
                                  output_directory = 'output',
                                  update_dimensions=True, extrabonds=True ):

1406
1407
1408
1409
1410
1411
1412
        num_steps = int(num_steps//12)*12
        minimization_steps = int(minimization_steps//24)*24
        if num_steps < 12:
            raise ValueError("Must run with at least 12  steps")
        if minimization_steps < 24:
            raise ValueError("Must run with at least 24 minimization steps")

1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
        format_data = self.__dict__.copy() # get parameters from System object

        format_data['extrabonds'] = """extraBonds on
extraBondsFile $prefix.exb
""" if extrabonds else ""

        if self.useTclForces:
            format_data['margin'] = ""
            format_data['tcl_forces'] = """tclForces on
tclForcesScript $prefix.forces.tcl
"""
        else:
            format_data['margin'] = """margin              30
"""
            format_data['tcl_forces'] = ""

        if update_dimensions:
1430
            format_data['dimensions'] = self.dimensions_from_structure()
1431

1432
        for k,v in zip('XYZ', format_data['dimensions']):
1433
1434
1435
            format_data['origin'+k] = -v*0.5
            format_data['cell'+k] = v

1436
        format_data['prefix'] = output_name
1437
1438
        format_data['minimization_steps'] = int(minimization_steps//2)
        format_data['num_steps'] = num_steps
1439
        format_data['output_directory'] = output_directory
1440
        filename = '{}.namd'.format(output_name)
1441
1442
1443
1444
1445
1446
1447

        with open(filename,'w') as fh:
            fh.write("""
set prefix {prefix}
set nLast 0;			# increment when continueing a simulation
set n [expr $nLast+1]
set out {output_directory}/$prefix-$n
1448
set temperature {temperature}
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520

structure          $prefix.psf
coordinates        $prefix.pdb

outputName         $out
XSTfile            $out.xst
DCDfile            $out.dcd

#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input
paraTypeCharmm	    on
parameters          charmm36.nbfix/par_all36_na.prm
parameters	    charmm36.nbfix/par_water_ions_na.prm

wrapAll             off

# Force-Field Parameters
exclude             scaled1-4
1-4scaling          1.0
switching           on
switchdist           8
cutoff              10
pairlistdist        12
{margin}

# Integrator Parameters
timestep            2.0  ;# 2fs/step
rigidBonds          all  ;# needed for 2fs steps
nonbondedFreq       1
fullElectFrequency  3
stepspercycle       12

# PME (for full-system periodic electrostatics)
PME                 no
PMEGridSpacing      1.2

# Constant Temperature Control
langevin            on    ;# do langevin dynamics
# langevinDamping     1   ;# damping coefficient (gamma); used in original study
langevinDamping     0.1   ;# less friction for faster relaxation
langevinTemp        $temperature
langevinHydrogen    off    ;# don't couple langevin bath to hydrogens

# output
useGroupPressure    yes
xstFreq             4800
outputEnergies      4800
dcdfreq             4800
restartfreq         48000

#############################################################
## EXTRA FORCES                                            ##
#############################################################

# ENM and intrahelical extrabonds
{extrabonds}
{tcl_forces}

#############################################################
## RUN                                                     ##
#############################################################

# Continuing a job from the restart files
cellBasisVector1 {cellX} 0 0
cellBasisVector2 0 {cellY} 0
cellBasisVector3 0 0 {cellZ}

if {{$nLast == 0}} {{
    temperature 300
1521
1522
1523
1524
    fixedAtoms on
    fixedAtomsForces on
    fixedAtomsFile $prefix.fixed.pdb
    fixedAtomsCol B
1525
    minimize {minimization_steps:d}
1526
    fixedAtoms off
1527
    minimize {minimization_steps:d}
1528
1529
1530
1531
1532
1533
1534
1535
}} else {{
    bincoordinates  {output_directory}/$prefix-$nLast.restart.coor
    binvelocities   {output_directory}/$prefix-$nLast.restart.vel
}}

run {num_steps:d}
""".format(**format_data))

1536
1537
    def atomic_simulate(self, output_name, output_directory='output', dry_run = False, namd2=None, log_file=None, num_procs=None, gpu=None, minimization_steps=4800, num_steps=1e6, write_pqr=False, copy_ff_from=get_resource_path("charmm36.nbfix") ):

cmaffeo2's avatar
cmaffeo2 committed
1538
1539
1540
1541
        if self.cacheUpToDate == False: # TODO: remove cache?
            self._countParticleTypes()
            self._updateParticleOrder()

1542
1543
        if output_directory == '': output_directory='.'
        self.writePdb( output_name + ".pdb" )
1544
        self.writePdb( output_name + ".fixed.pdb", beta_from_fixed=True )
cmaffeo2's avatar
cmaffeo2 committed
1545
        if write_pqr: self.write_pqr( output_name + ".pqr" )        
1546
        self.writePsf( output_name + ".psf" )
1547
        self.write_namd_configuration( output_name, output_directory = output_directory, minimization_steps=minimization_steps, num_steps=num_steps )
1548
1549
1550
1551
1552
1553
1554
1555

        if copy_ff_from is not None and copy_ff_from != '':
            try:
                shutil.copytree( copy_ff_from, Path(copy_ff_from).stem )
            except FileExistsError:
                pass

        
1556
        # os.sync()
1557

1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
        if not dry_run:
            if namd2 is None:
                for path in os.environ["PATH"].split(os.pathsep):
                    path = path.strip('"')
                    fname = os.path.join(path, "namd2")
                    if os.path.isfile(fname) and os.access(fname, os.X_OK):
                        namd2 = fname
                        break

            if namd2 is None: raise Exception("NAMD2 was not found")

            if not os.path.exists(namd2):
                raise Exception("NAMD2 was not found")
            if not os.path.isfile(namd2):
                raise Exception("NAMD2 was not found")
            if not os.access(namd2, os.X_OK):
                raise Exception("NAMD2 is not executable")

            if not os.path.exists(output_directory):
                os.makedirs(output_directory)
            elif not os.path.isdir(output_directory):
                raise Exception("output_directory '%s' is not a directory!" % output_directory)

            if num_procs is None:
                import multiprocessing
                num_procs = max(1,multiprocessing.cpu_count()-1)

            cmd = [namd2, '+p{}'.format(num_procs), "%s.namd" % output_name]
            cmd = tuple(str(x) for x in cmd)

            print("Running NAMD2 with: %s" % " ".join(cmd))
            if log_file is None or (hasattr(log_file,'write') and callable(log_file.write)):
                fd = sys.stdout if log_file is None else log_file
                process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
                for line in process.stdout:
                    fd.write(line)
                    fd.flush()
            else:
                with open(log_file,'w') as fd:
                    process = subprocess.Popen(cmd, stdout=log_file, universal_newlines=True)
                    process.communicate()
For faster browsing, not all history is shown. View entire blame