"vscode:/vscode.git/clone" did not exist on "b886644412b64749d3f184bfc5fcfc58c4e83d6c"
BrookVerletDynamics.cpp 26.6 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2008 Stanford University and the Authors.           *
 * Authors: Mark Friedrichs                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include <sstream>
#include "BrookVerletDynamics.h"
#include "BrookPlatform.h"
35
#include "openmm/OpenMMException.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
36
#include "BrookStreamImpl.h"
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
37
38
39
#include "kernels/kshakeh.h"
#include "kernels/kupdatemd.h"
#include "kernels/kcommon.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

using namespace OpenMM;
using namespace std;

/** 
 *
 * Constructor
 * 
 */

BrookVerletDynamics::BrookVerletDynamics( ){

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookVerletDynamics::BrookVerletDynamics";

56
57
58
   BrookOpenMMFloat zero                    = static_cast<BrookOpenMMFloat>(  0.0 );
   BrookOpenMMFloat one                     = static_cast<BrookOpenMMFloat>(  1.0 );
   BrookOpenMMFloat oneMinus                = static_cast<BrookOpenMMFloat>( -1.0 );
Mark Friedrichs's avatar
Mark Friedrichs committed
59
60
61

// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
62
   _numberOfParticles                 = -1;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
63
   _internalStepCount                 = 0;
Mark Friedrichs's avatar
Mark Friedrichs committed
64
65
66

   // mark stream dimension variables as unset

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
67
68
69
   _verletParticleStreamWidth         = -1;
   _verletParticleStreamHeight        = -1;
   _verletParticleStreamSize          = -1;
Mark Friedrichs's avatar
Mark Friedrichs committed
70
71
72
73
74

   for( int ii = 0; ii < LastStreamIndex; ii++ ){
      _verletStreams[ii]   = NULL;
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
75
   _stepSize                          = oneMinus;
Mark Friedrichs's avatar
Mark Friedrichs committed
76
77
78

   // setup inverse sqrt masses

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
79
   _inverseMasses                     = NULL;
Mark Friedrichs's avatar
Mark Friedrichs committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

}   
 
/** 
 * Destructor
 * 
 */

BrookVerletDynamics::~BrookVerletDynamics( ){

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookVerletDynamics::~BrookVerletDynamics";

// ---------------------------------------------------------------------------------------

   for( int ii = 0; ii < LastStreamIndex; ii++ ){
      delete _verletStreams[ii];
   }

   delete[] _inverseMasses;

}

/** 
 * Get stepSize
 * 
 * @return  stepSize
 *
 */

BrookOpenMMFloat BrookVerletDynamics::getStepSize( void ) const {
   return _stepSize;
}

/** 
 * Set stepSize
 * 
 * @param   stepSize
 *
 * @return      DefaultReturnValue
 *
 */

int BrookVerletDynamics::_setStepSize( BrookOpenMMFloat stepSize ){
   _stepSize = stepSize;
   return DefaultReturnValue;
}

/** 
 * Update parameters -- only way parameters can be set
 * 
 * @param  step size       step size
 *
 * @return   DefaultReturnValue
 *
 */

int BrookVerletDynamics::updateParameters( double stepSize ){

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
142
143
   static int showUpdate               = 1;
   static int maxShowUpdate            = 3;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
144
   static const std::string methodName = "BrookVerletDynamics::updateParameters";
Mark Friedrichs's avatar
Mark Friedrichs committed
145
146
147

   // ---------------------------------------------------------------------------------------

148
   _setStepSize( (BrookOpenMMFloat) stepSize );
Mark Friedrichs's avatar
Mark Friedrichs committed
149
150
151
152
153
154
   _updateVerletStreams( );

   // show update

   if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){
      std::string contents = getContentsString( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
155
      (void) fprintf( getLog(), "%s contents\n%s", methodName.c_str(), contents.c_str() );
Mark Friedrichs's avatar
Mark Friedrichs committed
156
157
158
159
160
161
      (void) fflush( getLog() );

   }

   return DefaultReturnValue;

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
162
}
Mark Friedrichs's avatar
Mark Friedrichs committed
163
164

/** 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
165
 * Get Particle stream size
Mark Friedrichs's avatar
Mark Friedrichs committed
166
 *
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
167
 * @return  Particle stream size
Mark Friedrichs's avatar
Mark Friedrichs committed
168
169
170
 *
 */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
171
172
int BrookVerletDynamics::getVerletDynamicsParticleStreamSize( void ) const {
   return _verletParticleStreamSize;
Mark Friedrichs's avatar
Mark Friedrichs committed
173
174
175
}

/** 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
176
 * Get particle stream width
Mark Friedrichs's avatar
Mark Friedrichs committed
177
 *
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
178
 * @return  particle stream width
Mark Friedrichs's avatar
Mark Friedrichs committed
179
180
181
 *
 */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
182
183
int BrookVerletDynamics::getVerletDynamicsParticleStreamWidth( void ) const {
   return _verletParticleStreamWidth;
Mark Friedrichs's avatar
Mark Friedrichs committed
184
185
186
}

/** 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
187
 * Get particle stream height
Mark Friedrichs's avatar
Mark Friedrichs committed
188
 *
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
189
 * @return particle stream height
Mark Friedrichs's avatar
Mark Friedrichs committed
190
191
 */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
192
193
int BrookVerletDynamics::getVerletDynamicsParticleStreamHeight( void ) const {
   return _verletParticleStreamHeight;
Mark Friedrichs's avatar
Mark Friedrichs committed
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
}

/** 
 * Get VPrime stream 
 *
 * @return  Vprime stream
 *
 */

BrookFloatStreamInternal* BrookVerletDynamics::getVPrimeStream( void ) const {
   return _verletStreams[VPrimeStream];
}

/** 
 * Get XPrime stream 
 *
 * @return  Xprime stream
 *
 */

BrookFloatStreamInternal* BrookVerletDynamics::getXPrimeStream( void ) const {
   return _verletStreams[XPrimeStream];
}

/** 
 * Get InverseMass stream 
 *
 * @return  inverse mass stream
 *
 */

BrookFloatStreamInternal* BrookVerletDynamics::getInverseMassStream( void ) const {
   return _verletStreams[InverseMassStream];
}

/** 
 * Initialize stream dimensions
 * 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
232
 * @param numberOfParticles             number of particles
Mark Friedrichs's avatar
Mark Friedrichs committed
233
234
235
236
237
238
 * @param platform                  platform
 *      
 * @return ErrorReturnValue if error, else DefaultReturnValue
 *
 */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
239
int BrookVerletDynamics::_initializeStreamSizes( int numberOfParticles, const Platform& platform ){
Mark Friedrichs's avatar
Mark Friedrichs committed
240
241
242
243
244
245
246

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookVerletDynamics::_initializeStreamSizes";

// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
247
248
249
   _verletParticleStreamSize     = getParticleStreamSize( platform );
   _verletParticleStreamWidth    = getParticleStreamWidth( platform );
   _verletParticleStreamHeight   = getParticleStreamHeight( platform );
Mark Friedrichs's avatar
Mark Friedrichs committed
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268

   return DefaultReturnValue;
}

/** 
 * Initialize streams
 * 
 * @param platform                  platform
 *
 * @return ErrorReturnValue if error, else DefaultReturnValue
 *
 */

int BrookVerletDynamics::_initializeStreams( const Platform& platform ){

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookVerletDynamics::_initializeStreams";

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
269
   BrookOpenMMFloat dangleValue         = (BrookOpenMMFloat) 0.0;
Mark Friedrichs's avatar
Mark Friedrichs committed
270
271
272

// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
273
274
   int sdParticleStreamSize             = getVerletDynamicsParticleStreamSize();
   int sdParticleStreamWidth            = getVerletDynamicsParticleStreamWidth();
Mark Friedrichs's avatar
Mark Friedrichs committed
275
276

    _verletStreams[VPrimeStream]        = new BrookFloatStreamInternal( BrookCommon::VPrimeStream,
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
277
278
                                                                        sdParticleStreamSize, sdParticleStreamWidth,
                                                                        BrookStreamInternal::Float3, dangleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
279
280

    _verletStreams[XPrimeStream]        = new BrookFloatStreamInternal( BrookCommon::XPrimeStream,
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
281
282
                                                                        sdParticleStreamSize, sdParticleStreamWidth,
                                                                        BrookStreamInternal::Float3, dangleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
283
284

    _verletStreams[InverseMassStream]   = new BrookFloatStreamInternal( BrookCommon::InverseMassStream,
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
285
286
                                                                        sdParticleStreamSize, sdParticleStreamWidth,
                                                                        BrookStreamInternal::Float, dangleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305

   return DefaultReturnValue;
}

/** 
 * Update sd streams -- called after parameters change
 * 
 * @return ErrorReturnValue if error, else DefaultReturnValue
 *
 */

int BrookVerletDynamics::_updateVerletStreams( void ){

// ---------------------------------------------------------------------------------------

   static const std::string methodName       = "BrookVerletDynamics::_updateVerletStreams";

// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
306
   int particleStreamSize                    = getVerletDynamicsParticleStreamSize();
Mark Friedrichs's avatar
Mark Friedrichs committed
307

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
308
309
   BrookOpenMMFloat* inverseMass = new BrookOpenMMFloat[particleStreamSize];
   memset( inverseMass, 0, particleStreamSize*sizeof( BrookOpenMMFloat ) ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
310

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
311
312
   int numberOfParticles                         = getNumberOfParticles();
   for( int ii = 0; ii < numberOfParticles; ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
313
314
315
316
317
318
319
320
321
322
323
324
325
326
      inverseMass[ii]  = _inverseMasses[ii];
   }

   _verletStreams[InverseMassStream]->loadFromArray( inverseMass );

   delete[] inverseMass;

   return DefaultReturnValue;

}

/** 
 * Set masses 
 * 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
327
 * @param masses             particle masses
Mark Friedrichs's avatar
Mark Friedrichs committed
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
 *
 */

int BrookVerletDynamics::_setInverseMasses( const std::vector<double>& masses ){

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookVerletDynamics::_setInverseSqrtMasses";

   BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
   BrookOpenMMFloat one                     = (BrookOpenMMFloat)  1.0;

// ---------------------------------------------------------------------------------------

   // setup inverse sqrt masses

   _inverseMasses = new BrookOpenMMFloat[masses.size()];
   int index          = 0;
   for( std::vector<double>::const_iterator ii = masses.begin(); ii != masses.end(); ii++, index++ ){
      if( *ii != 0.0 ){
         BrookOpenMMFloat value    = static_cast<BrookOpenMMFloat>(*ii);
         _inverseMasses[index] = ( one/value );
      } else {
         _inverseMasses[index] = zero;
      }
   }

   return DefaultReturnValue;
}   
 
/*  
 * Setup of VerletDynamics parameters
 *
 * @param masses                masses
 * @param platform              Brook platform
 *
 * @return nonzero value if error
 *
 * */
    
int BrookVerletDynamics::setup( const std::vector<double>& masses, const Platform& platform ){
    
// ---------------------------------------------------------------------------------------

   static const std::string methodName      = "BrookVerletDynamics::setup";

// ---------------------------------------------------------------------------------------

376
   const BrookPlatform& brookPlatform       = dynamic_cast<const BrookPlatform&> (platform);
Mark Friedrichs's avatar
Mark Friedrichs committed
377
378
   setLog( brookPlatform.getLog() );

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
379
380
   int numberOfParticles  = (int) masses.size();
   setNumberOfParticles( numberOfParticles );
Mark Friedrichs's avatar
Mark Friedrichs committed
381
382
383

   // set stream sizes and then create streams

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
384
   _initializeStreamSizes( numberOfParticles, platform );
Mark Friedrichs's avatar
Mark Friedrichs committed
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
   _initializeStreams( platform );

   _setInverseMasses( masses );

   return DefaultReturnValue;
}

/* 
 * Get contents of object
 *
 * @param level   level of dump
 *
 * @return string containing contents
 *
 * */

std::string BrookVerletDynamics::getContentsString( int level ) const {

// ---------------------------------------------------------------------------------------

   static const std::string methodName      = "BrookVerletDynamics::getContentsString";

   static const unsigned int MAX_LINE_CHARS = 256;
   char value[MAX_LINE_CHARS];
   static const char* Set                   = "Set";
   static const char* NotSet                = "Not set";

// ---------------------------------------------------------------------------------------

   std::stringstream message;
   std::string tab   = "   ";

#ifdef WIN32
#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
#else
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
#endif

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
423
424
   (void) LOCAL_SPRINTF( value, "%d", getNumberOfParticles() );
   message << _getLine( tab, "Number of particles:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
425

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
426
427
   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamWidth() );
   message << _getLine( tab, "Particle stream width:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
428

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
429
430
   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamHeight() );
   message << _getLine( tab, "Particle stream height:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
431

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
432
433
   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamSize() );
   message << _getLine( tab, "Particle stream size:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450

   (void) LOCAL_SPRINTF( value, "%.5f", getStepSize() );
   message << _getLine( tab, "Step size:", value ); 

   message << _getLine( tab, "Log:",                  (getLog()                ? Set : NotSet) ); 

   for( int ii = 0; ii < LastStreamIndex; ii++ ){
      message << std::endl;
      if( _verletStreams[ii] ){
         message << _verletStreams[ii]->getContentsString( );
      }
   }

#undef LOCAL_SPRINTF

   return message.str();
}
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
451

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
/** 
 * Update
 * 
 * @param  positions                   particle positions
 * @param  velocities                  particle velocities
 * @param  forces                      particle forces
 * @param  brookShakeAlgorithm         BrookShakeAlgorithm reference
 *
 * @return  DefaultReturnValue
 *
 */

int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImpl& velocityStream,
                                 const BrookStreamImpl& forceStreamC,
                                 BrookShakeAlgorithm& brookShakeAlgorithm ){

// ---------------------------------------------------------------------------------------

   static std::string methodName  = "\nBrookVerletDynamics::update";
Mark Friedrichs's avatar
Mark Friedrichs committed
471
   static int printOn             = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
472
   FILE* log;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
473
474
475

// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
476
477
   _internalStepCount++;

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
478
479
480
481
482
483
//setLog( stderr );
   if( printOn && getLog() ){
       log    = getLog();; 
   } else {
      printOn = 0;
   }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
484

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
485
486
   BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
487
   if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
488
489
490
491
492
493
494

      static int showAux = 1;
      if( showAux ){
         showAux = 0; 

/*
         std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
495
         (void) fprintf( log, "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
496
*/
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
497
498
         (void) fprintf( log, "%s step=%d \n%s\n\nShake contents\n%s", methodName.c_str(), _internalStepCount, getContentsString().c_str(),
                         brookShakeAlgorithm.getContentsString().c_str() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
499
         (void) fflush( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
      }    

   }

   // To Shake or not to Shake

   if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){

      // integration step
   
      kupdate_md1( (float) getStepSize(),
                   positionStream.getBrookStream(),
                   velocityStream.getBrookStream(),
                   forceStream.getBrookStream(),
                   getInverseMassStream()->getBrookStream(),
                   getXPrimeStream()->getBrookStream() 
                 );
    
      // diagnostics
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
520
521
522
      if( printOn ){
         (void) fprintf( log, "\n%s step=%d Post kupdate_md1: particleStrW=%3d step=%.5f",
                                   methodName.c_str(), _internalStepCount, getVerletDynamicsParticleStreamWidth(), getStepSize() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
523
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
524
525
         (void) fprintf( log, "\nInverseMassStream %d\n", _internalStepCount );
         getInverseMassStream()->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
526
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
527
         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
528
529
         (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
         brookStreamInternalPos->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
530
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
531
         (void) fprintf( log, "\nForceStream %d\n", _internalStepCount );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
532
         BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
533
534
535
536
537
538
539
540
541
         std::vector<std::vector<double> > forceStatistics;
         brookStreamInternalF->getStatistics( forceStatistics, getNumberOfParticles() );

         std::stringstream tag;
         tag << _internalStepCount << " Fxx ";
         std::string stats = brookStreamInternalF->printStatistics( tag.str(), forceStatistics );
         (void) fprintf( log, "\nStep %d Force stats:\n%s", _internalStepCount, stats.c_str() );

         brookStreamInternalF->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
542
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
543
         BrookStreamInternal* brookStreamInternalV   = velocityStream.getBrookStreamInternal();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
544
545
546
547
548
549
550
551
         std::vector<std::vector<double> > velocityStatistics;
         brookStreamInternalV->getStatistics( velocityStatistics, getNumberOfParticles() );
         std::stringstream tagV;
         tagV << _internalStepCount << " Vxx ";
         stats = brookStreamInternalPos->printStatistics( tagV.str(), velocityStatistics );
         (void) fprintf( log, "\nStep %d Velocity stats:\n%s", _internalStepCount, stats.c_str() );
         (void) fprintf( log, "\nVelocityStream %d\n", _internalStepCount );
         brookStreamInternalV->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
552
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
553
554
         (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
         getXPrimeStream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
      }   

      // Shake

      kshakeh_fix1( 
                    (float) brookShakeAlgorithm.getMaxIterations(),
                    (float) getVerletDynamicsParticleStreamWidth(),
                    brookShakeAlgorithm.getShakeTolerance(),
                    brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                    positionStream.getBrookStream(),
                    getXPrimeStream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
572
      if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
573

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
574
         (void) fprintf( log, "\n%s Post kshakeh_fix1: particleStrW=%3d", methodName.c_str(), getVerletDynamicsParticleStreamWidth() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
575
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
576
577
         (void) fprintf( log, "\nShakeParticleIndicesStream %d\n", _internalStepCount );
         brookShakeAlgorithm.getShakeParticleIndicesStream()->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
578
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
579
         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
580
581
         (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
         brookStreamInternalPos->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
582
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
583
584
         (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
         getXPrimeStream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
585

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
586
587
         (void) fprintf( log, "\nShakeParticleParameterStream %d\n", _internalStepCount );
         brookShakeAlgorithm.getShakeParticleParameterStream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
588

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
589
590
         (void) fprintf( log, "\nShakeXCons0\n" );
         brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
591

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
592
593
         (void) fprintf( log, "\nShakeXCons1\n" );
         brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
594

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
595
596
         (void) fprintf( log, "\nShakeXCons2\n" );
         brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
597

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
598
599
         (void) fprintf( log, "\nShakeXCons3\n" );
         brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
600
601
602
603
604
605

      }   

      // Shake gather
   
      kshakeh_update1_fix1( 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
606
                    (float) brookShakeAlgorithm.getShakeConstraintStreamWidth(),
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
607
608
609
610
611
612
613
614
615
                    brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
                    getXPrimeStream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
                    brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
                    getXPrimeStream()->getBrookStream() );
                    //positionStream.getBrookStream() );
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
616
      if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
617

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
618
619
         (void) fprintf( log, "\n%s Post kshakeh_update1_fix1: step=%d ShakeConstraintStreamWidth=%3d",
                                   methodName.c_str(), _internalStepCount, brookShakeAlgorithm.getShakeConstraintStreamWidth() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
620
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
621
622
         (void) fprintf( log, "\nShakeInverseMapStream %d\n", _internalStepCount );
         brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
623
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
624
         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
625
626
         (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
         brookStreamInternalPos->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
627
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
628
629
         (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
         getXPrimeStream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
630

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
631
632
         (void) fprintf( log, "\nShakeXCons0\n" );
         brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
633

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
634
635
         (void) fprintf( log, "\nShakeXCons1\n" );
         brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
636

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
637
638
         (void) fprintf( log, "\nShakeXCons2\n" );
         brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
639

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
640
641
         (void) fprintf( log, "\nShakeXCons3\n" );
         brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
642
643
644
645
646
647
648
649
650
651
652
653

      }   

      // second integration step

      float inverseStepSize = 1.0f/getStepSize();
      kupdate_md2( inverseStepSize,
                   getXPrimeStream()->getBrookStream(), 
                   positionStream.getBrookStream(),
                   velocityStream.getBrookStream(),
                   positionStream.getBrookStream() );

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
654
      if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
655

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
656
657
         (void) fprintf( log, "\n%s step=%d Post kupdate_md2: inverseStepSize=%3e",
                                   methodName.c_str(), _internalStepCount, inverseStepSize );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
658
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
659
         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
660
661
662
         std::string violationString;
         brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, violationString, 0.0001f );
         (void) fprintf( log, "Shake: %s\n", violationString.c_str() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
663
664
         (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
         brookStreamInternalPos->printToFile( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
665

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
666
667
         (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
         getXPrimeStream()->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
668

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
669
         brookStreamInternalPos = velocityStream.getBrookStreamInternal();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
670
671
         (void) fprintf( log, "\nVelocityStream %d\n", _internalStepCount );
         brookStreamInternalPos->printToFile( log ); 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
672
673
674
      }   

   } else {
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
675

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
676
677
678
679
680
681
682
683
684
      kupdateMdNoShake( getStepSize(),
                        positionStream.getBrookStream(),
                        velocityStream.getBrookStream(),
                        forceStream.getBrookStream(),
                        getInverseMassStream()->getBrookStream(),
                        velocityStream.getBrookStream(),
                        positionStream.getBrookStream() );
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
685
686
   // diagnostics

Mark Friedrichs's avatar
Mark Friedrichs committed
687
   if( 0 && (_internalStepCount % 100) == 0 ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
      FILE*	log1     = stderr;
      float  epsilon = 1.0e-01f;

      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
      BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamInternal();
      BrookStreamInternal* brookStreamInternalFrc  = forceStream.getBrookStreamInternal();

      // check for nan and infinities

      int coordinateNans                           = brookStreamInternalPos->checkForNans( );
      int velocityNans                             = brookStreamInternalVel->checkForNans( );
      int forceNans                                = brookStreamInternalFrc->checkForNans( );
      int abort                                    = abs( coordinateNans ) + abs( velocityNans ) + abs( forceNans );

      // Shake violations

      std::string violationString;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
705
706
707
708
709
710
      if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
         int constraintViolations                     = brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, violationString, 0.0001f );
         abort                                       += abs( constraintViolations );
      } else {
         violationString                              = "Shake not active";
      }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748

      // force sums ~ 0?
      std::vector<float> sums;
      brookStreamInternalFrc->sumColumns( sums );

      // check if should abort

      (void) fprintf( log1, "%d Nans: x=%d v=%d f=%d ", _internalStepCount, coordinateNans, velocityNans, forceNans );
      (void) fprintf( log1, " Fsum[" );
      for( int ii = 0; ii < 3; ii++ ){
         if( fabsf( sums[ii] ) > epsilon ){
            abort++;
         }
         (void) fprintf( log1, "%12.4e ", sums[ii] );
      }
      (void) fprintf( log1, "] %s abort=%d\n", violationString.c_str(), abort );
      (void) fflush( log1 );
      if( abort ){

         (void) fprintf( log1, "Aborting:\n" );

         brookStreamInternalPos->printToFile( log1 );
         brookStreamInternalVel->printToFile( log1 );
         brookStreamInternalFrc->printToFile( log1 );

         exit(1);
      }
  
/*
      std::vector<std::vector<double> > velocityStatistics;
      brookStreamInternalPos->getStatistics( velocityStatistics, getNumberOfParticles() );
      std::stringstream tagV;
      tagV << _internalStepCount << " Vxx "; 
      std::string stats = brookStreamInternalPos->printStatistics( tagV.str(), velocityStatistics );
      (void) fprintf( log1, "\nStep %d Velocity stats:\n%s", _internalStepCount, stats.c_str() );
*/
      //removeCom( brookStreamInternalPos, getInverseMassStream() );
   }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
749
750
751
752
753
//_brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );

   return DefaultReturnValue;

}