-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathIDPmd13_2.f90
2498 lines (1962 loc) · 81.2 KB
/
IDPmd13_2.f90
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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
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
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
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
376
377
378
379
380
381
382
383
384
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
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
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
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! *********************************************************
! molecular dynamics simulation on protein binding and folding
!
! author: Hueseyin Kaya revised by Zhirong Liu
!
! by Yongqi Huang
!
! SchemeO (0): chain A is moves without any influence from chain B (by ignoring interchain Q)
! SchemeA (1): chain B is fixed to avoid unfolding, but chain A is allowed to move freely
! SchemeB (2): both chains A and B are moving
!
! Periodic boundary condiction is used to make sure chain A is near chain B.
!
! ! Solvation model: M. S. Cheung et al, PNAS 99, 685 (2002)
! note that there is a little mistake for term C in the formula
!
! Simulate the system to see if it first reach the native
! or the denatural states
!
! restore when the data is running out
! ouput some conformation for further usage
!
!
! scale the bending, torsion, non-bonding interactions
! Non-native HP interactions with changable strength.
! Note that we use the appNCS file in data format of: residue_i, residue_j, contact_type (1: native contact; 2: non-native hydrophobic interaction; else: repulsive)
! A gap between residue number of two chains can be specified.
! Other columns than the first three are ignored.
! The order and number of contact intems can be random.
! Add the bias and histogram for Qw, which is defined as Qw=Qb-\alpha_w*\sum_native [r_AB-r_AB^(0)]
! The order and number of contact intems can be random.
! References:
! [1] Yongqi Huang and Zhirong Liu.
! Kinetic advantage of intrinsically disordered proteins in coupled folding-binding process: a critical assessment of the ¡°fly-casting¡± mechanism.
! J. Mol. Biol. 393 (5), 1143-1159 (2009).
! [2] Yongqi Huang and Zhirong Liu.
! Smoothing molecular interactions: the ¡°kinetic buffer¡± effect of intrinsically disordered proteins.
! Proteins: Struct. Funct. Bioinform. 78 (16), 3251-3259 (2010).
! [3] Yongqi Huang and Zhirong Liu.
! Nonnative interactions in coupled folding and binding processes of intrinsically disordered proteins.
! PLoS ONE 5 (11), e15375/1-10 (2010).
! [4] Huseyin Kaya and Hue Sun Chan.
! Solvation Effects and Driving Forces for Protein Thermodynamic and Kinetic Cooperativity: How Adequate is Native-centric Topological Modeling?
! J. Mol. Biol. 326 (3), 911-931 (2003).
! [5] Zhirong Liu and Hue Sun Chan.
! Solvation and desolvation effects in protein folding: native flexibility, kinetic cooperativity and enthalpic barriers under isostability conditions.
! Phys. Biol. 2 (4), S75-S85 (2005).
! **********************************************************
PROGRAM MOLECULARDYNAMICS
! *********************************************************
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
character*50 title, date, filename1, filename2, filename
character*50 term, errmess
parameter(MAXN=430,MAXUNBO=(MAXN)*(MAXN-1)/2)
parameter(nbinmax=105)
parameter(nEbinmax=105)
parameter(nRbinmax=105)
dimension xinit(MAXN),yinit(MAXN),zinit(MAXN)
dimension xsave(MAXN),ysave(MAXN),zsave(MAXN)
! number of particles for chain1, moving chains, chain1&2, and number of unbonding term (i<j)
! iFlagMov=0 move chain 1 without influece from chain 2; =1: move chain 1 only; =2: move both chains
common/Bnum/npart1, iFlagMov, nResGap, npartM, nparttol, nunbond
common/coordinates/x(MAXN),y(MAXN),z(MAXN)
common/velocity/vx(MAXN),vy(MAXN),vz(MAXN)
common/forcetotalnew/fx(MAXN),fy(MAXN),fz(MAXN)
common/forcetotalold/fxo(MAXN),fyo(MAXN),fzo(MAXN)
common/forcestretching/fxr(MAXN),fyr(MAXN),fzr(MAXN)
common/forcebending/fxth(MAXN),fyth(MAXN),fzth(MAXN)
common/forcetorsion/fxph(MAXN),fyph(MAXN),fzph(MAXN)
common/forceunbonded/fxun(MAXN),fyun(MAXN),fzun(MAXN)
common/forcerandnew/frandx(MAXN),frandy(MAXN),frandz(MAXN)
common/forcerandold/frandxo(MAXN),frandyo(MAXN),frandzo(MAXN)
common/nativeinfo/rbond_nat(MAXN),runbond_nat(MAXUNBO),theta_nat(MAXN),dihedral_nat(MAXN)
common/gyrationradius/gyr
common/constants/pi,boltz,avsn0,gamma,amass,sigma_ij
common/variables/temp,dt,nstep,nadim,nsnap,gm,enscale
common/contactmap/iun(MAXUNBO),jun(MAXUNBO),kunbond(MAXUNBO),nQnative_f1,nQnative_f2,nQnative_b
common/randomterms/c_0,c_1,c_2,randconst
common/interacparam/ck_r,ck_tht,ck_phi1,ck_phi3,epsil1,epsil2,epsil
common/LagrangeData/ga1_f1, ga2_f1, ga1_f2, ga2_f2, ga1_b,ga2_b,gQ0_f1,&
gQ0_f2, gQ0_b, gr0, gdr, gQ_f, gQ_f1, gQ_f2, gQ_b, eGr
common/LagrangeData2/ga1_w, ga2_w, gQ0_w, gQ_w, alpha_Qw
common/OutputConform/nConformOuput, outQf1_i, outQf1_f, outQf2_i, outQf2_f, iQb_i, iQb_f, nOutput0, ndOutput
! output nConfrom conformations with outQf1_i<=Qf1<=outQf1_f and outQf2_i<=Qf2<=outQf2_f and iQb_i<=Qb<=iQb_b
! noutput0, ndOutput: the initial step and gap step of the output
common/TransCoefSet/nConform,nCon0, nRunConf, gQbnativ, gQbdenatural, gQf1nativ, gQf1denatural, gQf2nativ, gQf2denatural
common/Setdtgm/ s_dt, s_gm, iFixKr
common/Solvation/k_sol, n_sol, m_sol, epsilon_p, epsilon_pp, pseudoQ_f,pseudoQ_b, dr_sol, ddr_sol
common/WithWithoutSolvation/ Is_Solvation
dimension ioctsd(4)
! histogram of (Q_f,Q_b), (Q_f), (Q_b)
common/HistogramData_fb/nbin_f,nbin_b,dbin_f,dbin_b,vbin0,nbinsnap,nbinsnap0
common/HistogramData_FB/PFBbin(nbinmax,nbinmax), PFbin(nbinmax), PBbin(nbinmax)
! histogram of (Q_f,Q_b,E)
common/EQ_fbHistogramData/nEbin,dEbin,vEbin0, IsEbin
common/EQ_FBHistogramData/PFBEbin(nbinmax,nbinmax,nEbinmax)
! histogram of (Q_b,E_b)
common/EbQbHistogramParameter/nEbbin,dEbbin,vEbbin0, IsEbbin
common/EbQbHistogramData/PQbEbbin(nbinmax,nEbinmax)
! histogram of (Q_f,Q_b,R)
common/RQ_fbHistogramData/nRbin,dRbin,vRbin0, IsRbin
common/RQ_FBHistogramData/PFBRbin(nbinmax,nbinmax,nRbinmax)
! histogram of (Q_f1, Q_f2, Q_b)
common/HistogramData_FFB/PFFBbin(nbinmax,nbinmax,nbinmax)
! histogram of (Q_w), (Q_w, Q_b) for construction of F(Q_b), (Q_w, R, Qb=0?) for F(R) and K, (Q_w, E_b, Qb=0?) for K changing with epsilon_b
common/Q_wHistogramParam/nWbin, dWbin, vWbin0, IsWbin, cri_Qb
common/Q_wHistogramData/PQwbin(nbinmax), PQwQbbin(nbinmax,nbinmax), PQwQfbin(nbinmax,nbinmax), &
PQwRbin(nbinmax,nRbinmax,2), PQwEbbin(nbinmax,nEbinmax,2)
! for periodic boundary condiction,R is the distance to the center of the box
common/pbc/pL,dl,R
! Rescaling the interaction strength of interchain interaction and intrachain interaction
common/bias/Alpha1,Alpha2,Beta,Delta
common/non_native/CritR_non, gQ_non_f, gQ_non_b
write(*, '('' Enter input file name: '')')
read (*, *) filename
open ( 100, file = filename, status = 'old' )
read ( 100, '( 2a14 )' ) title, date
! write( *, '( 2a14 )' ) title, date
! write( *, '()' )
read ( 100, *, end = 10 ) filename1
!-------------read chains -------------------
read ( 100, * ) npart1, nparttol, nResGap, iFlagMov
if(iFlagMov.eq.1 .or. iFlagMov.eq.0) then
npartM=npart1
else
npartM=nparttol
endif
if(npart1.gt.nparttol) call myErr('Err: npart1.gt.nparttol')
if(npart1.lt.1 .or. npart1.gt.MAXN) call myErr('Err: npart1.lt.1 .or. .gt.MAXN')
if(nparttol.gt.MAXN) call myErr('Err: nparttol.gt.MAXN')
! file of the initial conformation: (InitialConf.dat)
read ( 100, * ) filename2
open ( 20, file = filename2, status = 'old' )
! Native conformation: (NativeConf.dat)
read ( 100, * ) filename2
open ( 2, file = filename2, status = 'old' )
! contact map file: (appNCS.dat)
read ( 100, * ) filename2
open ( 4, file = filename2, status = 'old' )
read ( 100, * ) term
!-------------read parameters-------------------
read (100,'(4i4)') ioctsd
read (100,*) epsil, epsil1, epsil2, enscale
read (100,*) ck_r, ck_tht, ck_phi1, ck_phi3
read (100,*) sigma_ij, amass, gamma
read (100,*) avsn0, boltz
read (100,*) gr0, gdr
read (100,*) s_dt, s_gm, iFixKr
read (100,*) k_sol, n_sol, m_sol
read (100,*) epsilon_p, epsilon_pp
read (100,*) temp
read ( 100, * ) ga1_f1, ga2_f1, gQ0_f1
read ( 100, * ) ga1_f2, ga2_f2, gQ0_f2
read ( 100, * ) ga1_b, ga2_b, gQ0_b
read ( 100, * ) ga1_w, ga2_w, gQ0_w, alpha_Qw
read (100,*) Is_Solvation
read (100,*) nConform, nCon0, nRunConf
read (100,*) gQbnativ, gQbdenatural, gQf1nativ, gQf1denatural, gQf2nativ, gQf2denatural
read (100,*) nsnap, nstep
! read (100,*) ndOutput
read (100, *) nConformOutput, nOutput0, ndOutput
read (100, *) outQf1_i, outQf1_f, outQf2_i, outQf2_f, outQb_i, outQb_f
read (100,*) nbinsnap, nbinsnap0
read (100,*) nbin_f, nbin_b
read (100,*) dbin_f, dbin_b, vbin0
read (100,*) IsEbin, nEbin, dEbin, vEbin0
read (100,*) IsEbbin, nEbbin, dEbbin, vEbbin0
read (100,*) IsRbin, nRbin, dRbin, vRbin0
read (100,*) IsWbin, nWbin, dWbin, vWbin0, cri_Qb
read (100,*) pL, dl
read (100,*) Alpha1, Alpha2, Beta
read (100,*) Delta, CritR_non
if(nbin_f.gt.nbinmax) call myErr('Err: nbin_f.gt.nbinmax')
if(nbin_b.gt.nbinmax) call myErr('Err: nbin_b.gt.nbinmax')
if(nEbin.gt.nEbinmax) call myErr('Err: nEbin.gt.nEbinmax')
if(nRbin.gt.nRbinmax) call myErr('Err: nRbin.gt.nRbinmax')
!-------------open files to save-------------------
LF1 = len_trim(filename1)
! file to save output information: (Output*)
LF2 = LF1 + 6
filename2(1:6 ) = 'Output'
filename2(7:LF2) = filename1(1:LF1)
open ( 9, file = filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_b): (QbHist*)
LF2 = LF1 + 6
filename2(1:6 ) = 'QbHist'
filename2(7:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0) open (14,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_f): (QfHist*)
LF2 = LF1 + 6
filename2(1:6 ) = 'QfHist'
filename2(7:LF2) = filename1(1:LF1)
open (15,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_f,Q_b): (QfbHist*)
LF2 = LF1 + 7
filename2(1:7 ) = 'QfbHist'
filename2(8:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0) open (19,file=filename2(1:LF2), status = 'replace' )
! file to save the (Q,E)-histogram: (QEHist*.21)
LF2 = LF1 + 6
filename2(1:6 ) = 'QEHist'
filename2(7:LF2) = filename1(1:LF1)
if(IsEbin.eq.1) open (23,file=filename2(1:LF2), status = 'replace' )
! file to save the (Q,R)-histogram: (QRHist*.21)
LF2 = LF1 + 6
filename2(1:6 ) = 'QRHist'
filename2(7:LF2) = filename1(1:LF1)
if(IsRbin.eq.1) open (24,file=filename2(1:LF2), status = 'replace' )
! file to save the (Qf1,Qf2,Qb)-histogram: (QffbHist*.21)
LF2 = LF1 + 8
filename2(1:8 ) = 'QffbHist'
filename2(9:LF2) = filename1(1:LF1)
if(iFlagMov.eq.2) open (25,file=filename2(1:LF2), status = 'replace' )
! file to save the output conformation, d*.mol
LF2 = LF1 + 5
filename2(1:1) = 'd'
filename2(2 :LF1+1) = filename1(1:LF1)
filename2(LF1+2:LF2) = '.mol'
open (22,file=filename2(1:LF2), status = 'replace' )
! file to save the information for output conformations, d*.mol.dat
LF2 = LF1 + 9
filename2(1:1) = 'd'
filename2(2 :LF1+1) = filename1(1:LF1)
filename2(LF1+2:LF2) = '.mol.dat'
open(27,file=filename2(1:LF2), status = 'replace' )
! file to save the output conformation coordinate for further simulation, d*.X
LF2 = LF1 + 3
filename2(2 :LF1+1) = filename1(1:LF1)
filename2(1:1) = 'd'
filename2(LF1+2:LF2) = '.X'
open (26,file=filename2(1:LF2), status = 'replace' )
! file to save the binding/unbinding time and transmission coefficient, d*.tc
LF2 = LF1 + 4
filename2(2 :LF1+1) = filename1(1:LF1)
filename2(1:1) = 'd'
filename2(LF1+2:LF2) = '.tc'
open (28,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_b,E_b): (QbEbHist*)
LF2 = LF1 + 8
filename2(1:8 ) = 'QbEbHist'
filename2(9:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0 .and. IsEbbin.eq.1) open (29,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_w): (QwHist*)
LF2 = LF1 + 6
filename2(1:6 ) = 'QwHist'
filename2(7:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0 .and. IsWbin.eq.1) open (30,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_w,Q_b): (QwQbHist*)
LF2 = LF1 + 8
filename2(1:8 ) = 'QwQbHist'
filename2(9:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0 .and. IsWbin.eq.1) open (31,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_w,R): (QwRHist*)
LF2 = LF1 + 7
filename2(1:7 ) = 'QwRHist'
filename2(8:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0 .and. IsWbin.eq.1) open (32,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_w,E_b): (QwEbHist*)
LF2 = LF1 + 8
filename2(1:8 ) = 'QwEbHist'
filename2(9:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0 .and. IsWbin.eq.1) open (33,file=filename2(1:LF2), status = 'replace' )
! file to save Histogram P(Q_w,Q_f1): (QwEbHist*)
LF2 = LF1 + 8
filename2(1:8 ) = 'QwQfHist'
filename2(9:LF2) = filename1(1:LF1)
if(iFlagMov.gt.0 .and. IsWbin.eq.1) open (34,file=filename2(1:LF2), status = 'replace' )
write (*,'('' Length of chain 1 (npart1): '',t50, i10)') npart1
write (*,'('' Number of moving chains (iFlagMov): '',t50, i10)') iFlagMov
write (*,'('' Total length (nparttol): '',t50, i10)') nparttol
write (*,'('' Chain-Residue number gap (nResGap): '',t50, i10)') nResGap
write (*,'('' interaction strength epsilon'',t50, f10.5)') epsil
write (*,'('' interaction strength epsilon_1'',t50, f10.5)') epsil1
write (*,'('' interaction strength epsilon_2'',t50, f10.5)') epsil2
write (*,'('' interaction strength ck_r'',t50, f10.5)') ck_r
write (*,'('' interaction strength ck_tht'',t50, f10.5)') ck_tht
write (*,'('' interaction strength ck_phi1'',t50, f10.5)') ck_phi1
write (*,'('' interaction strength ck_phi3'',t50, f10.5)') ck_phi3
write (*,'('' interaction strength scaling param. '',t50, f10.5)')enscale
write (*,'('' Solvation k: '',t50, i10)') k_sol
write (*,'('' Solvation n: '',t50, i10)') n_sol
write (*,'('' Solvation m: '',t50, i10)') m_sol
write (*,'('' Solvation epsilon pri: '',t50, f10.5)') epsilon_p
write (*,'('' Solvation epsilon pripri:'',t50, f10.5)')epsilon_pp
write (*,'('' dt:'',t50,f10.6)') s_dt
write (*,'('' gamma:'',t50,f10.6)') s_gm
write (*,'('' Is Fix Kr:'',t50,i10)') iFixKr
write (*,'('' temperature'',t50, f10.5)') temp
write (*,'('' cut-off'',t50, f10.5)') gamma
write (*,'('' hard core distance'',t50, f10.5)') sigma_ij
write (*,'('' total simulation step'',t50, i13)') nstep
write (*,'('' frequency of snapshots '',t50, i10)') nsnap
write (*,'('' Avogadro number'',t50, f10.5)') avsn0
write (*,'('' Boltzmann constant'',t50, f10.5)') boltz
write (*,'('' random seed'',t50, 4i4)') ioctsd
write (*,'('' Lagrange constraint a1_f1: '',t50, f10.5)') ga1_f1
write (*,'('' Lagrange constraint a2_f1: '',t50, f10.5)') ga2_f1
write (*,'('' Lagrange constraint Q0_f1: '',t50, f10.5)') gQ0_f1
write (*,'('' Lagrange constraint a1_f2: '',t50, f10.5)') ga1_f2
write (*,'('' Lagrange constraint a2_f2: '',t50, f10.5)') ga2_f2
write (*,'('' Lagrange constraint Q0_f2: '',t50, f10.5)') gQ0_f2
write (*,'('' Lagrange constraint a1_b: '',t50, f10.5)') ga1_b
write (*,'('' Lagrange constraint a2_b: '',t50, f10.5)') ga2_b
write (*,'('' Lagrange constraint Q0_b: '',t50, f10.5)') gQ0_b
write (*,'('' Lagrange constraint a1_b: '',t50, f10.5)') ga1_2
write (*,'('' Lagrange constraint a2_b: '',t50, f10.5)') ga2_2
write (*,'('' Lagrange constraint Q0_b: '',t50, f10.5)') gQ0_2
write (*,'('' Parameter for Q_w as alpha_Qw: '',t50, f10.5)') alpha_Qw
write (*,'('' Lagrange constraint r0: '',t50, f10.5)') gr0
write (*,'('' Lagrange constraint dr: '',t50, f10.5)') gdr
write (*,'('' nConformOutput:'',t50,i10)') nConformOutput
write (*,'('' nOutput0 :'',t50,i10)') nOutput0
write (*,'('' ndOutput:'',t50,i10)') ndOutput
write (*,'('' outQf1_i: '',t50, f10.5)') outQf1_i
write (*,'('' outQf1_f: '',t50, f10.5)') outQf1_f
write (*,'('' outQf2_i: '',t50, f10.5)') outQf2_i
write (*,'('' outQf2_f: '',t50, f10.5)') outQf2_f
write (*,'('' outQb_i: '',t50, f10.5)') outQb_i
write (*,'('' outQb_f: '',t50, f10.5)') outQb_f
write (*,'('' Histogram nbin_f: '',t50, i10)') nbin_f
write (*,'('' Histogram dbin_f: '',t50, f10.5)') dbin_f
write (*,'('' Histogram nbin_b: '',t50, i10)') nbin_b
write (*,'('' Histogram dbin_b: '',t50, f10.5)') dbin_b
write (*,'('' num. of conformation:'',t50,i10)')nConform
write (*,'('' num. of conformation to skip:'',t50,i10)')nCon0
write (*,'('' nRunConf'',t50,i10)') nRunConf
write (*,'('' gQbnativ'',t50,f10.3)') gQbnativ
write (*,'('' gQbdenatural:'',t50,f10.3)') gQbdenatural
write (*,'('' gQf1nativ'',t50,f10.3)') gQf1nativ
write (*,'('' gQf1denatural:'',t50,f10.3)') gQf1denatural
write (*,'('' gQf2nativ'',t50,f10.3)') gQf2nativ
write (*,'('' gQf2denatural:'',t50,f10.3)') gQf2denatural
write (*,'('' Force scheme: '',t50,i10)') Is_Solvation
write (*,'('' If save E-Histogram: '',t50,i10)') IsEbin
write (*,'('' E-Histogram nbin: '',t50, i10)') nEbin
write (*,'('' E-Histogram dbin: '',t50, f10.5)') dEbin
write (*,'('' E-Histogram vbin0: '',t50, f10.5)') vEbin0
write (*,'('' If save R-Histogram: '',t50,i10)') IsRbin
write (*,'('' R-Histogram nbin: '',t50, i10)') nRbin
write (*,'('' R-Histogram dbin: '',t50, f10.5)') dRbin
write (*,'('' R-Histogram vbin0: '',t50, f10.5)') vRbin0
write (*,'('' If save Qw-Histogram: '',t50,i10)') IsWbin
write (*,'('' Qw-Histogram nbin: '',t50, i10)') nWbin
write (*,'('' Qw-Histogram dbin: '',t50, f10.5)') dWbin
write (*,'('' Qw-Histogram vbin0: '',t50, f10.5)') vWbin0
write (*,'('' cri_Qb for Qw-Histogram: '',t50, f10.5)') cri_Qb
write (*,'('' periodic boundary condiction size:'',t50, f10.5)') pL
write (*,'('' periodic boundary condiction buffer:'',t50, f10.5)') dl
write (*,'('' Ratio of intrachain interaction 1:'',t50, f10.5)') Alpha1
write (*,'('' Ratio of intrachain interaction 2:'',t50, f10.5)') Alpha2
write (*,'('' Ratio of interchain interaction:'',t50, f10.5)') Beta
write (*,'('' Delta:'',t50, f10.5)') Delta
write (*,'('' CritR_non:'',t50, f10.5)') CritR_non
read (2,*)(x(j),y(j),z(j),j=1,nparttol)
nunbond=0
do while(nunbond.lt.MAXUNBO)
read(4,*, ERR=17, END=17) item1, item2, item3
if(item1.gt.npart1) then
if(item1.lt.npart1+nResGap) call myErr('contact map error 4.')
item1=item1-nResGap
endif
if(item2.gt.npart1) then
if(item2.lt.npart1+nResGap) call myErr('contact map error 5.')
item2=item2-nResGap
endif
if(item1.le.0.or.item1.gt.nparttol) call myErr('contact map error 1.')
if(item2.le.0.or.item2.gt.nparttol) call myErr('contact map error 2.')
if(item1.eq.item2) call myErr('contact map error 3.')
nunbond=nunbond+1
if(item1.lt.item2) then
iun(nunbond)=item1
jun(nunbond)=item2
else
iun(nunbond)=item2
jun(nunbond)=item1
endif
kunbond(nunbond)=item3
if(item1.gt.npartM .and. item2.gt.npartM .or. &
iFlagMov.eq.0 .and. (item1.gt.npart1 .or. item2.gt.npart1)) nunbond=nunbond-1
end do
17 continue
write (*,'('' Effective Number of unbond: '',i7)') nunbond
! read (4,*)( iun(j), jun(j), kunbond(j), runbond_nat(j), j = 1, nunbond )
call setrn(ioctsd)
call nativeinformation
call intpar(enerkin)
write (*,'('' Native Qf:'',i4,i4,''; Qb:'',i4)') nQnative_f1,nQnative_f2,nQnative_b
! total number of folding/unfolding events
nFoldTol=0
nunFoldTol=0
aveNadim=0
do 316 j=1,nCon0
do 315 i=1,nparttol
read(20,*) xinit(i),yinit(i),zinit(i)
315 continue
316 continue
!===================== simulation start =====================
do 2000 nConfcount=nCon0+1, nConform
do 320 i=1,nparttol
read(20,*) xinit(i),yinit(i),zinit(i)
320 continue
nFoldSub=0
nunFoldSub=0
do 1500 nRuncount=1, nRunConf
do 330 i=1, nparttol
x(i)=xinit(i)
y(i)=yinit(i)
z(i)=zinit(i)
330 continue
call origin_adjust
call InitVel(enerkin)
call force(e_pot,e_unbond_tot,e_bind_tot,e_tors_tot,e_bend_tot,e_bond_tot)
call RANTERM
do 319 kl=1,npartM
fxo(kl)=fx(kl)
fyo(kl)=fy(kl)
fzo(kl)=fz(kl)
frandxo(kl)=frandx(kl)
frandyo(kl)=frandy(kl)
frandzo(kl)=frandz(kl)
319 continue
nOutputCount=0
nOutputGap=-1
nadim=-1
! main dynamics cycle
do 100 nadim_new=0,nstep
nadim=nadim+1
call verlet(enerkin,e_pot)
!-------------save data for purpose of restore-------------
if(mod(nadim,nsnap).eq.0 .and. vx(1).ge.-1e6 .and. vx(1).le.1e6) then
do 40 kl=1,nparttol
xsave(kl)=x(kl)
ysave(kl)=y(kl)
zsave(kl)=z(kl)
40 continue
nadim_old=nadim
nOutGap_old=nOutputGap
endif
! restore
if(.not.(vx(1).ge.-1e6 .and. vx(1).le.1e6)) then
do 45 i=1, nparttol
x(i)=xsave(i)
y(i)=ysave(i)
z(i)=zsave(i)
45 continue
call InitVel(enerkin)
call force(e_pot,e_unbond_tot,e_bind_tot,e_tors_tot,e_bend_tot,e_bond_tot)
call RANTERM
do 48 kl=1,npartM
fxo(kl)=fx(kl)
fyo(kl)=fy(kl)
fzo(kl)=fz(kl)
frandxo(kl)=frandx(kl)
frandyo(kl)=frandy(kl)
frandzo(kl)=frandz(kl)
48 continue
nadim=nadim_old
nOutputGap=nOutGap_old
endif
!-------------save data for purpose of restore end----------
!--------------output conformation-------------
nOutputGap=nOutputGap+1
if ( nadim == 0 ) write ( 27, '( '' nOutputCount, nadim, gQ_f, gQ_b, R'' )' )
if(nOutputCount.lt.nConformOutput .and. nadim.ge.nOutput0 .and. nOutputGap.ge.ndOutput .and. &
gQ_f1.ge.outQf1_i .and. gQ_f1.le.outQf1_f .and. &
gQ_f2.ge.outQf2_i .and. gQ_f2.le.outQf2_f .and. &
gQ_b.ge.outQb_i .and. gQ_b.le.outQb_f) then
nOutputGap=0
call write_mol(nOutputCount)
write ( 27, '( i7, i13, 3f10.3 )' ) nOutputCount, nadim, gQ_f, gQ_b, R
do 50 i=1, nparttol
write(26, *) x(i),y(i),z(i)
50 continue
write(26,'('' '')')
write(*,'(''n='',i3,'', nadim='',i8,'', gQ='',f7.3, '', old gQ='', 1i3)')nOutputCount,nadim,gQ,natcont
endif
!--------------output conformation end---------
if(gQ_f1.le.gQf1denatural .and. gQ_f2.le.gQf2denatural .and. gQ_b.le.gQbdenatural &
.or. gQ_f1.ge.gQf1nativ .and. gQ_f2.ge.gQf2nativ .and. gQ_b.ge.gQbnativ) then
goto 1000
endif
100 continue
1000 continue
aveNadim=aveNadim+nadim
if(.not.(vx(1).ge.-1e8 .and. vx(1).le.1e8)) then
k=0
else if(gQ_f1.le.gQf1denatural .and. gQ_f2.le.gQf2denatural .and. gQ_b.le.gQbdenatural) then
k=-1
nunFoldSub=nunFoldSub+1
nunFoldTol=nunFoldTol+1
else if(gQ_f1.ge.gQf1nativ .and. gQ_f2.ge.gQf2nativ .and. gQ_b.ge.gQbnativ) then
k=1
nFoldSub=nFoldSub+1
nFoldTol=nFoldTol+1
else
k=0
endif
write(28,'(3i5,i13,3f10.3)') nConfcount, nRuncount, k, nadim, gQ_b, gQ_f1, gQ_f2
1500 continue
write(28,'(''# Conform No.:'',i5,'', Folding times:'',i5,'',&
Transmission coefficient:'',f10.3,'', unFolding times:'',i5)') &
nConfcount,nFoldSub,1.0*nFoldSub/nRunConf, nunFoldSub
2000 continue
aveNadim=aveNadim/((nConform-nCon0)*nRunConf)
write(28,'(''# Total folding times:'',i5,'', Total unfolding times:'',i5)') nFoldTol, nunFoldTol
write(28,'(''# Average Running Time :'',f14.3)') aveNadim
!===================== simulation end =====================
! output target
call write_target(nOutputCount)
! output histogram
call write_histogram
close(2)
close(4)
close(20)
close(5)
close(9)
if(iFlagMov.gt.0) close(14)
close(15)
if(iFlagMov.gt.0) close(19)
close(21)
close(22)
if(IsEbin.eq.1) close(23)
if(IsRbin.eq.1) close(24)
if(iFlagMov.eq.2) close(25)
close(26)
close(27)
close(28)
if(iFlagMov.gt.0 .and. IsWbin.eq.1) then
close(30)
close(31)
close(32)
close(33)
close(34)
endif
10 continue
end
! *****************************************************
SUBROUTINE myErr(errmess)
! ****************************************************
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
character errmess(100)
write(*,*) errmess
write(9,*) errmess
call exit(0)
return
end
! *******************************************
SUBROUTINE write_histogram
! *******************************************
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
parameter(nbinmax=105)
parameter(nEbinmax=105)
parameter(nRbinmax=105)
common/Bnum/npart1, iFlagMov, nResGap, npartM, nparttol, nunbond
! histogram of (Q_f,Q_b)
common/HistogramData_fb/nbin_f,nbin_b,dbin_f,dbin_b,vbin0,nbinsnap,nbinsnap0
common/HistogramData_FB/PFBbin(nbinmax,nbinmax), PFbin(nbinmax), PBbin(nbinmax)
! histogram of (Q_f,Q_b,E)
common/EQ_fbHistogramData/nEbin,dEbin,vEbin0, IsEbin
common/EQ_FBHistogramData/PFBEbin(nbinmax,nbinmax,nEbinmax)
! histogram of (Q_b,E_b)
common/EbQbHistogramParameter/nEbbin,dEbbin,vEbbin0, IsEbbin
common/EbQbHistogramData/PQbEbbin(nbinmax,nEbinmax)
! histogram of (Q_f,Q_b,R)
common/RQ_fbHistogramData/nRbin,dRbin,vRbin0, IsRbin
common/RQ_FBHistogramData/PFBRbin(nbinmax,nbinmax,nRbinmax)
! histogram of (Q_f1, Q_f2, Q_b)
common/HistogramData_FFB/PFFBbin(nbinmax,nbinmax,nbinmax)
! histogram of (Q_w), (Q_w, Q_b) for construction of F(Q_b), (Q_w, R, Qb=0?) for F(R) and K, (Q_w, E_b, Qb=0?) for K changing with epsilon_b
common/Q_wHistogramParam/nWbin, dWbin, vWbin0, IsWbin, cri_Qb
common/Q_wHistogramData/PQwbin(nbinmax), PQwQbbin(nbinmax,nbinmax), PQwQfbin(nbinmax,nbinmax), &
PQwRbin(nbinmax,nRbinmax,2), PQwEbbin(nbinmax,nEbinmax,2)
! histogram result saving to file
do 20 i=1,nbin_f
write(15, '(f6.2, f12.1)') vbin0+(i-0.5)*dbin_f, PFbin(i)
20 continue
if (iFlagMov.le.0) goto 25
do 1001 i=1,nbin_b
sumQ_b=0
do 1002 j=1,nbin_f
sumQ_b=sumQ_b+PFBbin(j,i)
1002 continue
write(14, '(f6.2, f12.1)') vbin0+(i-0.5)*dbin_b,sumQ_b
1001 continue
25 continue
do 1003 i=1, nbin_f
do 1004 j=1, nbin_b
if(iFlagMov.gt.0 .and. PFBbin(i,j).gt.1e-5) write(19,'(2f7.2, f12.1)') vbin0 &
+(i-0.5)*dbin_f,vbin0+(j-0.5)*dbin_b,PFBbin(i,j)
if ( IsEbin == 0 ) goto 1007
do 1005 k=1,nEbin
if(PFBEbin(i,j,k).gt.1e-5) then
write(23, '(2f7.2, f8.2, f12.1)') vbin0+(i-0.5)*dbin_f,vbin0+(j-0.5)*dbin_b,vEbin0+(k-0.5)*dEbin,PFBEbin(i,j,k)
endif
1005 continue
1007 continue
if ( IsRbin == 0 ) goto 1008
do 1006 l=1,nRbin
if(PFBRbin(i,j,l).gt.1e-5) then
write(24, '(3f7.2, f12.1)') vbin0+(i-0.5)*dbin_f,vbin0+(j-0.5)*dbin_b,vRbin0+(l-0.5)*dRbin,PFBRbin(i,j,l)
endif
1006 continue
1008 continue
if ( iFlagMov .ne. 2 ) goto 1010
do 1009 l=1,nbin_f
if(PFFBbin(i,l,j).gt.1e-5) then
write(25, '(3f7.2, f12.1)') vbin0+(i-0.5)*dbin_f,vbin0+(l-0.5)*dbin_f,vbin0+(j-0.5)*dbin_b,PFFBbin(i,l,j)
endif
1009 continue
1010 continue
1004 continue
1003 continue
if ( iFlagMov.le.0 .or. IsEbbin.eq.0 ) goto 40
do 35 j=1, nbin_b
do 30 k=1,nEbbin
if(PQbEbbin(j,k).gt.1e-5) then
write(29, '(2f8.2, f12.1)') vbin0+(j-0.5)*dbin_b, vEbbin0+(k-0.5)*dEbbin,PQbEbbin(j,k)
endif
30 continue
35 continue
40 continue
if ( iFlagMov.le.0 .or. IsWbin.eq.0 ) goto 180
do 150 i=1, nWbin
write(30, '(f8.2, f12.1)') vWbin0+(i-0.5)*dWbin, PQwbin(i)
do 110 j=1, nbin_b
if(PQwQbbin(i,j).gt.1e-5) then
write(31, '(2f8.2, f12.1)') vWbin0+(i-0.5)*dWbin, vbin0+(j-0.5)*dbin_b, PQwQbbin(i,j)
endif
110 continue
do 120 l=1,nRbin
if(PQwRbin(i,l,1).gt.1e-5 .or. PQwRbin(i,l,2).gt.1e-5) then
write(32, '(2f8.2, 2f12.1)') vWbin0+(i-0.5)*dWbin, vRbin0+(l-0.5)*dRbin, PQwRbin(i,l,1), PQwRbin(i,l,2)
endif
120 continue
do 130 k=1,nEbbin
if(PQwEbbin(i,k,1).gt.1e-5 .or. PQwEbbin(i,k,2).gt.1e-5) then
write(33, '(2f8.2, 2f12.1)') vWbin0+(i-0.5)*dWbin, vEbbin0+(k-0.5)*dEbbin, PQwEbbin(i,k,1), PQwEbbin(i,k,2)
endif
130 continue
do 140 j=1, nbin_f
if(PQwQfbin(i,j).gt.1e-5) then
write(34, '(2f8.2, f12.1)') vWbin0+(i-0.5)*dWbin, vbin0+(j-0.5)*dbin_f, PQwQfbin(i,j)
endif
140 continue
150 continue
180 continue
end
! *******************************************
SUBROUTINE write_mol(nOutputCount)
! output conformations in MOL format
! *******************************************
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
parameter(MAXN=430,MAXUNBO=(MAXN)*(MAXN-1)/2)
common/Bnum/npart1, iFlagMov, nResGap, npartM, nparttol, nunbond
common/coordinates/x(MAXN),y(MAXN),z(MAXN)
nOutputCount = nOutputCount + 1
write(22,'(''Molecule-'',i7)') nOutputCount
write(22,'('' ViewerPro 3D 0'',/)')
if(iFlagMov.eq.1 .or. iFlagMov.eq.0) then
write(22,'(2i3,'' 0 0 0 0 0 0 0 0999 V2000'')') npart1,npart1-1
else
write(22,'(2i3,'' 0 0 0 0 0 0 0 0999 V2000'')') nparttol,nparttol-2
endif
do 50 i=1, npartM
write(22, '(3f10.4,'' C 0 0 0 0 0 0 0 0 0 0'')') x(i),y(i),z(i)
50 continue
do 211 i=1,npartM-1
j=i+1
if(j.le.npart1.or.i.gt.npart1) write(22,'(2i3,'' 1 0 0 0'')')i,j
211 continue
write(22,'(''M END'')')
write(22,'(''$$$$'')')
return
end
! *******************************************
SUBROUTINE write_target(nOutputCount)
! output conformations in MOL format
! *******************************************
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
parameter(MAXN=430,MAXUNBO=(MAXN)*(MAXN-1)/2)
common/Bnum/npart1, iFlagMov, nResGap, npartM, nparttol, nunbond
common/coordinates/x(MAXN),y(MAXN),z(MAXN)
nOutputCount=nOutputCount+1
write(22,'(''Molecule-'',i7)') nOutputCount
write(22,'('' ViewerPro 3D 0'',/)')
write(22,'(2i3,'' 0 0 0 0 0 0 0 0999 V2000'')') nparttol-npart1,nparttol-npart1-1
do 511 i=npart1+1, nparttol
write(22, '(3f10.4,'' C 0 0 0 0 0 0 0 0 0 0'')') x(i),y(i),z(i)
511 continue
do 512 i=1,nparttol-npart1-1
j=i+1
write(22,'(2i3,'' 1 0 0 0'')')i,j
512 continue
write(22,'(''M END'')')
write(22,'(''$$$$'')')
return
end
! *******************************************
SUBROUTINE origin_adjust
! place the target to (0,0,0)
! *******************************************
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
parameter(MAXN=430,MAXUNBO=(MAXN)*(MAXN-1)/2)
common/Bnum/npart1, iFlagMov, nResGap, npartM, nparttol, nunbond
common/coordinates/x(MAXN),y(MAXN),z(MAXN)
if(nparttol-npart1.le.0) goto 300
x00=0.0
y00=0.0
z00=0.0
do 281 j=npart1+1,nparttol
x00=x00+x(j)
y00=y00+y(j)
z00=z00+z(j)
281 continue
x00=x00/(nparttol-npart1)
y00=y00/(nparttol-npart1)
z00=z00/(nparttol-npart1)
do 282 j=1,nparttol
x(j)=x(j)-x00
y(j)=y(j)-y00
z(j)=z(j)-z00
282 continue
300 continue
end
! *******************************************
SUBROUTINE INTPAR(enerkin)
! *******************************************
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
parameter(MAXN=430,MAXUNBO=(MAXN)*(MAXN-1)/2)
common/Bnum/npart1, iFlagMov, nResGap, npartM, nparttol, nunbond
parameter(nbinmax=105)
parameter(nEbinmax=105)
parameter(nRbinmax=105)
common/velocity/vx(MAXN),vy(MAXN),vz(MAXN)
common/constants/pi,boltz,avsn0,gamma,amass,sigma_ij
common/variables/temp,dt,nstep,nadim,nsnap,gm,enscale
common/randomterms/c_0,c_1,c_2,randconst
common/interacparam/ck_r,ck_tht,ck_phi1,ck_phi3,epsil1,epsil2,epsil
common/Setdtgm/ s_dt, s_gm, iFixKr
common/Solvation/k_sol, n_sol, m_sol, epsilon_p, epsilon_pp, pseudoQ_f,pseudoQ_b, dr_sol, ddr_sol
! histogram of (Q_f,Q_b)
common/HistogramData_fb/nbin_f,nbin_b,dbin_f,dbin_b,vbin0,nbinsnap,nbinsnap0
common/HistogramData_FB/PFBbin(nbinmax,nbinmax), PFbin(nbinmax), PBbin(nbinmax)
! histogram of (Q_f,Q_b,E)
common/EQ_fbHistogramData/nEbin,dEbin,vEbin0, IsEbin
common/EQ_FBHistogramData/PFBEbin(nbinmax,nbinmax,nEbinmax)
! histogram of (Q_b,E_b)
common/EbQbHistogramParameter/nEbbin,dEbbin,vEbbin0, IsEbbin
common/EbQbHistogramData/PQbEbbin(nbinmax,nEbinmax)
! histogram of (Q_f,Q_b,R)
common/RQ_fbHistogramData/nRbin,dRbin,vRbin0, IsRbin
common/RQ_FBHistogramData/PFBRbin(nbinmax,nbinmax,nRbinmax)
! histogram of (Q_f1, Q_f2, Q_b)
common/HistogramData_FFB/PFFBbin(nbinmax,nbinmax,nbinmax)
! histogram of (Q_w), (Q_w, Q_b) for construction of F(Q_b), (Q_w, R, Qb=0?) for F(R) and K, (Q_w, E_b, Qb=0?) for K changing with epsilon_b
common/Q_wHistogramParam/nWbin, dWbin, vWbin0, IsWbin, cri_Qb
common/Q_wHistogramData/PQwbin(nbinmax), PQwQbbin(nbinmax,nbinmax), PQwQfbin(nbinmax,nbinmax), &
PQwRbin(nbinmax,nRbinmax,2), PQwEbbin(nbinmax,nEbinmax,2)
common/non_native/CritR_non, gQ_non_f, gQ_non_b
temp=temp*epsil/boltz
timeunit=sqrt(amass/epsil)*sigma_ij
! dt=0.005*timeunit
! gm=0.05/timeunit
dt=s_dt*timeunit
gm=s_gm/timeunit
c_0=dt*gm/(2.0*amass)
c_1=(1.0-c_0)*(1.0-c_0+c_0**2)
c_2=(dt/(2*amass))*(1.0-c_0+c_0**2)
randconst=sqrt(2.0*boltz*temp*gm/dt)
! no change ck_r -ZRLiu
! ck_r = enscale*ck_r
if(iFixKr.eq.0) ck_r=enscale*ck_r
ck_tht = enscale*ck_tht
ck_phi1= enscale*ck_phi1
ck_phi3= enscale*ck_phi3
epsil1 = enscale*epsil1
epsil2 = enscale*epsil2
epsilon_p=enscale*epsilon_p
epsilon_pp=enscale*epsilon_pp
! ccccccccccccccccccccccccccccccccc
! initiliaze velocities
! ccccccccccccccccccccccccccccccccc
sumvx=0.0
sumvy=0.0
sumvz=0.0
vm=sqrt(temp*boltz/amass)
do 1001 i=1,npartM
vx(i)=vm*gauss(xsi)
vy(i)=vm*gauss(xsi)
vz(i)=vm*gauss(xsi)
sumvx=sumvx+vx(i)
sumvy=sumvy+vy(i)
sumvz=sumvz+vz(i)
1001 continue