2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 call readi(controlcard,'SYM',symetr,1)
102 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
103 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
104 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
105 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
106 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
107 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
108 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
109 call reada(controlcard,'DRMS',drms,0.1D0)
110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
111 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
112 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
113 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
114 write (iout,'(a,f10.1)')'DRMS = ',drms
115 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
116 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
118 call readi(controlcard,'NZ_START',nz_start,0)
119 call readi(controlcard,'NZ_END',nz_end,0)
120 call readi(controlcard,'IZ_SC',iz_sc,0)
122 safety = 60.0d0*safety
125 call reada(controlcard,"T_BATH",t_bath,300.0d0)
126 minim=(index(controlcard,'MINIMIZE').gt.0)
127 dccart=(index(controlcard,'CART').gt.0)
128 overlapsc=(index(controlcard,'OVERLAP').gt.0)
129 overlapsc=.not.overlapsc
130 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
131 searchsc=.not.searchsc
132 sideadd=(index(controlcard,'SIDEADD').gt.0)
133 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
134 outpdb=(index(controlcard,'PDBOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 AFMlog=(index(controlcard,'AFM'))
141 print *,'AFMlog',AFMlog
142 call readi(controlcard,'IPRINT',iprint,0)
143 call readi(controlcard,'MAXGEN',maxgen,10000)
144 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
145 call readi(controlcard,"KDIAG",kdiag,0)
146 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
147 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
148 & write (iout,*) "RESCALE_MODE",rescale_mode
149 split_ene=index(controlcard,'SPLIT_ENE').gt.0
150 if (index(controlcard,'REGULAR').gt.0.0D0) then
151 call reada(controlcard,'WEIDIS',weidis,0.1D0)
155 if (index(controlcard,'CHECKGRAD').gt.0) then
157 if (index(controlcard,'CART').gt.0) then
159 elseif (index(controlcard,'CARINT').gt.0) then
164 elseif (index(controlcard,'THREAD').gt.0) then
166 call readi(controlcard,'THREAD',nthread,0)
167 if (nthread.gt.0) then
168 call reada(controlcard,'WEIDIS',weidis,0.1D0)
171 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
172 stop 'Error termination in Read_Control.'
174 else if (index(controlcard,'MCMA').gt.0) then
176 else if (index(controlcard,'MCEE').gt.0) then
178 else if (index(controlcard,'MULTCONF').gt.0) then
180 else if (index(controlcard,'MAP').gt.0) then
182 call readi(controlcard,'MAP',nmap,0)
183 else if (index(controlcard,'CSA').gt.0) then
185 crc else if (index(controlcard,'ZSCORE').gt.0) then
187 crc ZSCORE is rm from UNRES, modecalc=9 is available
190 cfcm else if (index(controlcard,'MCMF').gt.0) then
192 else if (index(controlcard,'SOFTREG').gt.0) then
194 else if (index(controlcard,'CHECK_BOND').gt.0) then
196 else if (index(controlcard,'TEST').gt.0) then
198 else if (index(controlcard,'MD').gt.0) then
200 else if (index(controlcard,'RE ').gt.0) then
204 lmuca=index(controlcard,'MUCA').gt.0
205 call readi(controlcard,'MUCADYN',mucadyn,0)
206 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
207 if (lmuca .and. (me.eq.king .or. .not.out1file ))
209 write (iout,*) 'MUCADYN=',mucadyn
210 write (iout,*) 'MUCASMOOTH=',muca_smooth
213 iscode=index(controlcard,'ONE_LETTER')
214 indphi=index(controlcard,'PHI')
215 indback=index(controlcard,'BACK')
216 iranconf=index(controlcard,'RAND_CONF')
217 i2ndstr=index(controlcard,'USE_SEC_PRED')
218 gradout=index(controlcard,'GRADOUT').gt.0
219 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
220 C DISTCHAINMAX become obsolete for periodic boundry condition
221 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
222 C Reading the dimensions of box in x,y,z coordinates
223 call reada(controlcard,'BOXX',boxxsize,100.0d0)
224 call reada(controlcard,'BOXY',boxysize,100.0d0)
225 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
226 c Cutoff range for interactions
227 call reada(controlcard,"R_CUT",r_cut,15.0d0)
228 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
229 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
230 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
231 if (lipthick.gt.0.0d0) then
232 bordliptop=(boxzsize+lipthick)/2.0
233 bordlipbot=bordliptop-lipthick
235 if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0))
236 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
237 buflipbot=bordlipbot+lipbufthick
238 bufliptop=bordliptop-lipbufthick
239 if ((lipbufthick*2.0d0).gt.lipthick)
240 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
242 write(iout,*) "bordliptop=",bordliptop
243 write(iout,*) "bordlipbot=",bordlipbot
244 write(iout,*) "bufliptop=",bufliptop
245 write(iout,*) "buflipbot=",buflipbot
248 if (me.eq.king .or. .not.out1file )
249 & write (iout,*) "DISTCHAINMAX",distchainmax
251 if(me.eq.king.or..not.out1file)
252 & write (iout,'(2a)') diagmeth(kdiag),
253 & ' routine used to diagonalize matrices.'
256 c--------------------------------------------------------------------------
257 subroutine read_REMDpar
261 implicit real*8 (a-h,o-z)
263 include 'COMMON.IOUNITS'
264 include 'COMMON.TIME1'
267 include 'COMMON.LANGEVIN'
269 include 'COMMON.LANGEVIN.lang0'
271 include 'COMMON.INTERACT'
272 include 'COMMON.NAMES'
274 include 'COMMON.REMD'
275 include 'COMMON.CONTROL'
276 include 'COMMON.SETUP'
278 character*320 controlcard
279 character*3200 controlcard1
280 integer iremd_m_total
282 if(me.eq.king.or..not.out1file)
283 & write (iout,*) "REMD setup"
285 call card_concat(controlcard)
286 call readi(controlcard,"NREP",nrep,3)
287 call readi(controlcard,"NSTEX",nstex,1000)
288 call reada(controlcard,"RETMIN",retmin,10.0d0)
289 call reada(controlcard,"RETMAX",retmax,1000.0d0)
290 mremdsync=(index(controlcard,'SYNC').gt.0)
291 call readi(controlcard,"NSYN",i_sync_step,100)
292 restart1file=(index(controlcard,'REST1FILE').gt.0)
293 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
294 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
295 if(max_cache_traj_use.gt.max_cache_traj)
296 & max_cache_traj_use=max_cache_traj
297 if(me.eq.king.or..not.out1file) then
298 cd if (traj1file) then
299 crc caching is in testing - NTWX is not ignored
300 cd write (iout,*) "NTWX value is ignored"
301 cd write (iout,*) " trajectory is stored to one file by master"
302 cd write (iout,*) " before exchange at NSTEX intervals"
304 write (iout,*) "NREP= ",nrep
305 write (iout,*) "NSTEX= ",nstex
306 write (iout,*) "SYNC= ",mremdsync
307 write (iout,*) "NSYN= ",i_sync_step
308 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
311 if (index(controlcard,'TLIST').gt.0) then
313 call card_concat(controlcard1)
314 read(controlcard1,*) (remd_t(i),i=1,nrep)
315 if(me.eq.king.or..not.out1file)
316 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
319 if (index(controlcard,'MLIST').gt.0) then
321 call card_concat(controlcard1)
322 read(controlcard1,*) (remd_m(i),i=1,nrep)
323 if(me.eq.king.or..not.out1file) then
324 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
327 iremd_m_total=iremd_m_total+remd_m(i)
329 write (iout,*) 'Total number of replicas ',iremd_m_total
332 if(me.eq.king.or..not.out1file)
333 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
336 c--------------------------------------------------------------------------
337 subroutine read_MDpar
341 implicit real*8 (a-h,o-z)
343 include 'COMMON.IOUNITS'
344 include 'COMMON.TIME1'
347 include 'COMMON.LANGEVIN'
349 include 'COMMON.LANGEVIN.lang0'
351 include 'COMMON.INTERACT'
352 include 'COMMON.NAMES'
354 include 'COMMON.SETUP'
355 include 'COMMON.CONTROL'
356 include 'COMMON.SPLITELE'
358 character*320 controlcard
360 call card_concat(controlcard)
361 call readi(controlcard,"NSTEP",n_timestep,1000000)
362 call readi(controlcard,"NTWE",ntwe,100)
363 call readi(controlcard,"NTWX",ntwx,1000)
364 call reada(controlcard,"DT",d_time,1.0d-1)
365 call reada(controlcard,"DVMAX",dvmax,2.0d1)
366 call reada(controlcard,"DAMAX",damax,1.0d1)
367 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
368 call readi(controlcard,"LANG",lang,0)
369 RESPA = index(controlcard,"RESPA") .gt. 0
370 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
371 ntime_split0=ntime_split
372 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
373 ntime_split0=ntime_split
374 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
375 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
376 rest = index(controlcard,"REST").gt.0
377 tbf = index(controlcard,"TBF").gt.0
378 usampl = index(controlcard,"USAMPL").gt.0
380 mdpdb = index(controlcard,"MDPDB").gt.0
381 call reada(controlcard,"T_BATH",t_bath,300.0d0)
382 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
383 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
384 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
385 if (count_reset_moment.eq.0) count_reset_moment=1000000000
386 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
387 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
388 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
389 if (count_reset_vel.eq.0) count_reset_vel=1000000000
390 large = index(controlcard,"LARGE").gt.0
391 print_compon = index(controlcard,"PRINT_COMPON").gt.0
392 rattle = index(controlcard,"RATTLE").gt.0
393 c if performing umbrella sampling, fragments constrained are read from the fragment file
399 if(me.eq.king.or..not.out1file) then
401 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
403 write (iout,'(a)') "The units are:"
404 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
405 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
406 & " acceleration: angstrom/(48.9 fs)**2"
407 write (iout,'(a)') "energy: kcal/mol, temperature: K"
409 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
410 write (iout,'(a60,f10.5,a)')
411 & "Initial time step of numerical integration:",d_time,
413 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
415 write (iout,'(2a,i4,a)')
416 & "A-MTS algorithm used; initial time step for fast-varying",
417 & " short-range forces split into",ntime_split," steps."
418 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
419 & r_cut," lambda",rlamb
421 write (iout,'(2a,f10.5)')
422 & "Maximum acceleration threshold to reduce the time step",
423 & "/increase split number:",damax
424 write (iout,'(2a,f10.5)')
425 & "Maximum predicted energy drift to reduce the timestep",
426 & "/increase split number:",edriftmax
427 write (iout,'(a60,f10.5)')
428 & "Maximum velocity threshold to reduce velocities:",dvmax
429 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
430 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
431 if (rattle) write (iout,'(a60)')
432 & "Rattle algorithm used to constrain the virtual bonds"
436 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
437 call reada(controlcard,"RWAT",rwat,1.4d0)
438 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
439 surfarea=index(controlcard,"SURFAREA").gt.0
440 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
441 if(me.eq.king.or..not.out1file)then
442 write (iout,'(/a,$)') "Langevin dynamics calculation"
445 & " with direct integration of Langevin equations"
446 else if (lang.eq.2) then
447 write (iout,'(a/)') " with TINKER stochasic MD integrator"
448 else if (lang.eq.3) then
449 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
450 else if (lang.eq.4) then
451 write (iout,'(a/)') " in overdamped mode"
453 write (iout,'(//a,i5)')
454 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
457 write (iout,'(a60,f10.5)') "Temperature:",t_bath
458 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
459 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
460 write (iout,'(a60,f10.5)')
461 & "Scaling factor of the friction forces:",scal_fric
462 if (surfarea) write (iout,'(2a,i10,a)')
463 & "Friction coefficients will be scaled by solvent-accessible",
464 & " surface area every",reset_fricmat," steps."
466 c Calculate friction coefficients and bounds of stochastic forces
467 eta=6*pi*cPoise*etawat
468 if(me.eq.king.or..not.out1file)
469 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
471 gamp=scal_fric*(pstok+rwat)*eta
472 stdfp=dsqrt(2*Rb*t_bath/d_time)
474 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
475 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
477 if(me.eq.king.or..not.out1file)then
478 write (iout,'(/2a/)')
479 & "Radii of site types and friction coefficients and std's of",
480 & " stochastic forces of fully exposed sites"
481 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
483 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
484 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
488 if(me.eq.king.or..not.out1file)then
489 write (iout,'(a)') "Berendsen bath calculation"
490 write (iout,'(a60,f10.5)') "Temperature:",t_bath
491 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
493 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
494 & count_reset_moment," steps"
496 & write (iout,'(a,i10,a)')
497 & "Velocities will be reset at random every",count_reset_vel,
501 if(me.eq.king.or..not.out1file)
502 & write (iout,'(a31)') "Microcanonical mode calculation"
504 if(me.eq.king.or..not.out1file)then
505 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
507 write(iout,*) "MD running with constraints."
508 write(iout,*) "Equilibration time ", eq_time, " mtus."
509 write(iout,*) "Constraining ", nfrag," fragments."
510 write(iout,*) "Length of each fragment, weight and q0:"
512 write (iout,*) "Set of restraints #",iset
514 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
515 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
517 write(iout,*) "constraints between ", npair, "fragments."
518 write(iout,*) "constraint pairs, weights and q0:"
520 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
521 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
523 write(iout,*) "angle constraints within ", nfrag_back,
524 & "backbone fragments."
525 write(iout,*) "fragment, weights:"
527 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
528 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
529 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
532 iset=mod(kolor,nset)+1
535 if(me.eq.king.or..not.out1file)
536 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
539 c------------------------------------------------------------------------------
542 C Read molecular data.
544 implicit real*8 (a-h,o-z)
550 include 'COMMON.IOUNITS'
553 include 'COMMON.INTERACT'
554 include 'COMMON.LOCAL'
555 include 'COMMON.NAMES'
556 include 'COMMON.CHAIN'
557 include 'COMMON.FFIELD'
558 include 'COMMON.SBRIDGE'
559 include 'COMMON.HEADER'
560 include 'COMMON.CONTROL'
561 include 'COMMON.DBASE'
562 include 'COMMON.THREAD'
563 include 'COMMON.CONTACTS'
564 include 'COMMON.TORCNSTR'
565 include 'COMMON.TIME1'
566 include 'COMMON.BOUNDS'
568 include 'COMMON.SETUP'
569 character*4 sequence(maxres)
571 double precision x(maxvar)
572 character*256 pdbfile
573 character*320 weightcard
574 character*80 weightcard_t,ucase
575 dimension itype_pdb(maxres)
576 common /pizda/ itype_pdb
577 logical seq_comp,fail
578 double precision energia(0:n_ene)
584 C Read weights of the subsequent energy terms.
585 call card_concat(weightcard)
586 call reada(weightcard,'WLONG',wlong,1.0D0)
587 call reada(weightcard,'WSC',wsc,wlong)
588 call reada(weightcard,'WSCP',wscp,wlong)
589 call reada(weightcard,'WELEC',welec,1.0D0)
590 call reada(weightcard,'WVDWPP',wvdwpp,welec)
591 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
592 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
593 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
594 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
595 call reada(weightcard,'WTURN3',wturn3,1.0D0)
596 call reada(weightcard,'WTURN4',wturn4,1.0D0)
597 call reada(weightcard,'WTURN6',wturn6,1.0D0)
598 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
599 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
600 call reada(weightcard,'WBOND',wbond,1.0D0)
601 call reada(weightcard,'WTOR',wtor,1.0D0)
602 call reada(weightcard,'WTORD',wtor_d,1.0D0)
603 call reada(weightcard,'WANG',wang,1.0D0)
604 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
605 call reada(weightcard,'SCAL14',scal14,0.4D0)
606 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
607 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
608 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
609 call reada(weightcard,'TEMP0',temp0,300.0d0)
610 call reada(weightcard,'WLT',wliptran,0.0D0)
611 if (index(weightcard,'SOFT').gt.0) ipot=6
612 C 12/1/95 Added weight for the multi-body term WCORR
613 call reada(weightcard,'WCORRH',wcorr,1.0D0)
614 if (wcorr4.gt.0.0d0) wcorr=wcorr4
634 if(me.eq.king.or..not.out1file)
635 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
636 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
638 10 format (/'Energy-term weights (unscaled):'//
639 & 'WSCC= ',f10.6,' (SC-SC)'/
640 & 'WSCP= ',f10.6,' (SC-p)'/
641 & 'WELEC= ',f10.6,' (p-p electr)'/
642 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
643 & 'WBOND= ',f10.6,' (stretching)'/
644 & 'WANG= ',f10.6,' (bending)'/
645 & 'WSCLOC= ',f10.6,' (SC local)'/
646 & 'WTOR= ',f10.6,' (torsional)'/
647 & 'WTORD= ',f10.6,' (double torsional)'/
648 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
649 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
650 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
651 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
652 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
653 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
654 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
655 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
656 & 'WTURN6= ',f10.6,' (turns, 6th order)')
657 if(me.eq.king.or..not.out1file)then
658 if (wcorr4.gt.0.0d0) then
659 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
660 & 'between contact pairs of peptide groups'
661 write (iout,'(2(a,f5.3/))')
662 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
663 & 'Range of quenching the correlation terms:',2*delt_corr
664 else if (wcorr.gt.0.0d0) then
665 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
666 & 'between contact pairs of peptide groups'
668 write (iout,'(a,f8.3)')
669 & 'Scaling factor of 1,4 SC-p interactions:',scal14
670 write (iout,'(a,f8.3)')
671 & 'General scaling factor of SC-p interactions:',scalscp
673 r0_corr=cutoff_corr-delt_corr
675 aad(i,1)=scalscp*aad(i,1)
676 aad(i,2)=scalscp*aad(i,2)
677 bad(i,1)=scalscp*bad(i,1)
678 bad(i,2)=scalscp*bad(i,2)
680 call rescale_weights(t_bath)
681 if(me.eq.king.or..not.out1file)
682 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
683 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
685 22 format (/'Energy-term weights (scaled):'//
686 & 'WSCC= ',f10.6,' (SC-SC)'/
687 & 'WSCP= ',f10.6,' (SC-p)'/
688 & 'WELEC= ',f10.6,' (p-p electr)'/
689 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
690 & 'WBOND= ',f10.6,' (stretching)'/
691 & 'WANG= ',f10.6,' (bending)'/
692 & 'WSCLOC= ',f10.6,' (SC local)'/
693 & 'WTOR= ',f10.6,' (torsional)'/
694 & 'WTORD= ',f10.6,' (double torsional)'/
695 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
696 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
697 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
698 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
699 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
700 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
701 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
702 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
703 & 'WTURN6= ',f10.6,' (turns, 6th order)')
704 if(me.eq.king.or..not.out1file)
705 & write (iout,*) "Reference temperature for weights calculation:",
707 call reada(weightcard,"D0CM",d0cm,3.78d0)
708 call reada(weightcard,"AKCM",akcm,15.1d0)
709 call reada(weightcard,"AKTH",akth,11.0d0)
710 call reada(weightcard,"AKCT",akct,12.0d0)
711 call reada(weightcard,"V1SS",v1ss,-1.08d0)
712 call reada(weightcard,"V2SS",v2ss,7.61d0)
713 call reada(weightcard,"V3SS",v3ss,13.7d0)
714 call reada(weightcard,"EBR",ebr,-5.50D0)
715 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
717 dyn_ss_mask(i)=.false.
721 dyn_ssbond_ij(i,j)=1.0d300
724 call reada(weightcard,"HT",Ht,0.0D0)
726 ss_depth=ebr/wsc-0.25*eps(1,1)
727 Ht=Ht/wsc-0.25*eps(1,1)
728 akcm=akcm*wstrain/wsc
729 akth=akth*wstrain/wsc
730 akct=akct*wstrain/wsc
731 v1ss=v1ss*wstrain/wsc
732 v2ss=v2ss*wstrain/wsc
733 v3ss=v3ss*wstrain/wsc
735 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
738 if(me.eq.king.or..not.out1file) then
739 write (iout,*) "Parameters of the SS-bond potential:"
740 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
742 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
743 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
744 write (iout,*)" HT",Ht
745 print *,'indpdb=',indpdb,' pdbref=',pdbref
747 if (indpdb.gt.0 .or. pdbref) then
748 read(inp,'(a)') pdbfile
749 if(me.eq.king.or..not.out1file)
750 & write (iout,'(2a)') 'PDB data will be read from file ',
751 & pdbfile(:ilen(pdbfile))
752 open(ipdbin,file=pdbfile,status='old',err=33)
754 33 write (iout,'(a)') 'Error opening PDB file.'
757 c print *,'Begin reading pdb data'
759 c print *,'Finished reading pdb data'
760 if(me.eq.king.or..not.out1file)
761 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
762 & ' nstart_sup=',nstart_sup
764 itype_pdb(i)=itype(i)
768 nct=nstart_sup+nsup-1
769 call contact(.false.,ncont_ref,icont_ref,co)
772 if(me.eq.king.or..not.out1file)
773 & write(iout,*)'Adding sidechains'
777 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
780 do while (fail.and.nsi.le.maxsi)
781 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
784 if(fail) write(iout,*)'Adding sidechain failed for res ',
785 & i,' after ',nsi,' trials'
790 if (indpdb.eq.0) then
791 C Read sequence if not taken from the pdb file.
793 c print *,'nres=',nres
794 if (iscode.gt.0) then
795 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
797 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
799 C Convert sequence to numeric code
801 itype(i)=rescode(i,sequence(i),iscode)
803 C Assign initial virtual bond lengths
809 vbld(i+nres)=dsc(iabs(itype(i)))
810 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
811 c write (iout,*) "i",i," itype",itype(i),
812 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
816 c print '(20i4)',(itype(i),i=1,nres)
819 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
821 if (itype(i).eq.ntyp1) then
825 else if (iabs(itype(i+1)).ne.20) then
827 else if (iabs(itype(i)).ne.20) then
834 if(me.eq.king.or..not.out1file)then
835 write (iout,*) "ITEL"
837 write (iout,*) i,itype(i),itel(i)
839 print *,'Call Read_Bridge.'
842 C 8/13/98 Set limits to generating the dihedral angles
847 read (inp,*) ndih_constr
848 if (ndih_constr.gt.0) then
850 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
851 if(me.eq.king.or..not.out1file)then
853 & 'There are',ndih_constr,' constraints on phi angles.'
855 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
859 phi0(i)=deg2rad*phi0(i)
860 drange(i)=deg2rad*drange(i)
862 if(me.eq.king.or..not.out1file)
863 & write (iout,*) 'FTORS',ftors
866 phibound(1,ii) = phi0(i)-drange(i)
867 phibound(2,ii) = phi0(i)+drange(i)
874 write (iout,'(a)') 'Boundaries in phi angle sampling:'
876 write (iout,'(a3,i5,2f10.1)')
877 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
883 cd print *,'NNT=',NNT,' NCT=',NCT
884 if (itype(1).eq.ntyp1) nnt=2
885 if (itype(nres).eq.ntyp1) nct=nct-1
887 if(me.eq.king.or..not.out1file)
888 & write (iout,'(a,i3)') 'nsup=',nsup
890 if (nsup.le.(nct-nnt+1)) then
891 do i=0,nct-nnt+1-nsup
892 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
898 & 'Error - sequences to be superposed do not match.'
901 do i=0,nsup-(nct-nnt+1)
902 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
904 nstart_sup=nstart_sup+i
910 & 'Error - sequences to be superposed do not match.'
913 if (nsup.eq.0) nsup=nct-nnt
914 if (nstart_sup.eq.0) nstart_sup=nnt
915 if (nstart_seq.eq.0) nstart_seq=nnt
916 if(me.eq.king.or..not.out1file)
917 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
918 & ' nstart_seq=',nstart_seq
920 c--- Zscore rms -------
921 if (nz_start.eq.0) nz_start=nnt
922 if (nz_end.eq.0 .and. nsup.gt.0) then
924 else if (nz_end.eq.0) then
927 if(me.eq.king.or..not.out1file)then
928 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
929 write (iout,*) 'IZ_SC=',iz_sc
931 c----------------------
934 if (.not.pdbref) then
935 call read_angles(inp,*38)
937 38 write (iout,'(a)') 'Error reading reference structure.'
939 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
940 stop 'Error reading reference structure'
944 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
954 call contact(.true.,ncont_ref,icont_ref,co)
956 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
958 if (constr_dist.gt.0) call read_dist_constr
959 write (iout,*) "After read_dist_constr nhpb",nhpb
960 if (AFMlog.gt.0) call read_afminp
962 if(me.eq.king.or..not.out1file)
963 & write (iout,*) 'Contact order:',co
965 if(me.eq.king.or..not.out1file)
966 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
969 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
971 if(me.eq.king.or..not.out1file)
972 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
973 & icont_ref(1,i),' ',
974 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
978 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
979 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
980 & modecalc.ne.10) then
981 C If input structure hasn't been supplied from the PDB file read or generate
983 if (iranconf.eq.0 .and. .not. extconf) then
984 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
985 & write (iout,'(a)') 'Initial geometry will be read in.'
987 read(inp,'(8f10.5)',end=36,err=36)
988 & ((c(l,k),l=1,3),k=1,nres),
989 & ((c(l,k+nres),l=1,3),k=nnt,nct)
990 write (iout,*) "Exit READ_CART"
991 write (iout,'(8f10.5)')
992 & ((c(l,k),l=1,3),k=1,nres),
993 & ((c(l,k+nres),l=1,3),k=nnt,nct)
994 call int_from_cart1(.true.)
995 write (iout,*) "Finish INT_TO_CART"
998 dc(j,i)=c(j,i+1)-c(j,i)
999 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1003 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1005 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1006 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1012 call read_angles(inp,*36)
1015 36 write (iout,'(a)') 'Error reading angle file.'
1017 call mpi_finalize( MPI_COMM_WORLD,IERR )
1019 stop 'Error reading angle file.'
1021 else if (extconf) then
1022 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1023 & write (iout,'(a)') 'Extended chain initial geometry.'
1025 theta(i)=90d0*deg2rad
1028 phi(i)=180d0*deg2rad
1031 alph(i)=110d0*deg2rad
1034 omeg(i)=-120d0*deg2rad
1035 if (itype(i).le.0) omeg(i)=-omeg(i)
1038 if(me.eq.king.or..not.out1file)
1039 & write (iout,'(a)') 'Random-generated initial geometry.'
1043 if (me.eq.king .or. fg_rank.eq.0 .and. (
1044 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1048 call gen_rand_conf(itmp,*30)
1050 30 write (iout,*) 'Failed to generate random conformation',
1051 & ', itrial=',itrial
1052 write (*,*) 'Processor:',me,
1053 & ' Failed to generate random conformation',
1063 write (iout,'(a,i3,a)') 'Processor:',me,
1064 & ' error in generating random conformation.'
1065 write (*,'(a,i3,a)') 'Processor:',me,
1066 & ' error in generating random conformation.'
1069 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1074 & ' error in generating random conformation.'
1079 elseif (modecalc.eq.4) then
1080 read (inp,'(a)') intinname
1081 open (intin,file=intinname,status='old',err=333)
1082 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1083 & write (iout,'(a)') 'intinname',intinname
1084 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1086 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1088 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1090 stop 'Error opening angle file.'
1094 C Generate distance constraints, if the PDB structure is to be regularized.
1095 if (nthread.gt.0) then
1096 call read_threadbase
1099 if (me.eq.king .or. .not. out1file)
1101 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1102 write (iout,'(/a,i3,a)')
1103 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1104 write (iout,'(20i4)') (iss(i),i=1,ns)
1106 write(iout,*)"Running with dynamic disulfide-bond formation"
1108 write (iout,'(/a/)') 'Pre-formed links are:'
1114 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1115 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1121 if (ns.gt.0.and.dyn_ss) then
1125 forcon(i-nss)=forcon(i)
1132 dyn_ss_mask(iss(i))=.true.
1135 if (i2ndstr.gt.0) call secstrp2dihc
1136 c call geom_to_var(nvar,x)
1137 c call etotal(energia(0))
1138 c call enerprint(energia(0))
1139 c call briefout(0,etot)
1141 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1142 cd write (iout,'(a)') 'Variable list:'
1143 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1145 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1146 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1147 & 'Processor',myrank,': end reading molecular data.'
1151 c--------------------------------------------------------------------------
1152 logical function seq_comp(itypea,itypeb,length)
1154 integer length,itypea(length),itypeb(length)
1157 if (itypea(i).ne.itypeb(i)) then
1165 c-----------------------------------------------------------------------------
1166 subroutine read_bridge
1167 C Read information about disulfide bridges.
1168 implicit real*8 (a-h,o-z)
1169 include 'DIMENSIONS'
1173 include 'COMMON.IOUNITS'
1174 include 'COMMON.GEO'
1175 include 'COMMON.VAR'
1176 include 'COMMON.INTERACT'
1177 include 'COMMON.LOCAL'
1178 include 'COMMON.NAMES'
1179 include 'COMMON.CHAIN'
1180 include 'COMMON.FFIELD'
1181 include 'COMMON.SBRIDGE'
1182 include 'COMMON.HEADER'
1183 include 'COMMON.CONTROL'
1184 include 'COMMON.DBASE'
1185 include 'COMMON.THREAD'
1186 include 'COMMON.TIME1'
1187 include 'COMMON.SETUP'
1188 C Read bridging residues.
1189 read (inp,*) ns,(iss(i),i=1,ns)
1191 if(me.eq.king.or..not.out1file)
1192 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1193 C Check whether the specified bridging residues are cystines.
1195 if (itype(iss(i)).ne.1) then
1196 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1197 & 'Do you REALLY think that the residue ',
1198 & restyp(itype(iss(i))),i,
1199 & ' can form a disulfide bridge?!!!'
1200 write (*,'(2a,i3,a)')
1201 & 'Do you REALLY think that the residue ',
1202 & restyp(itype(iss(i))),i,
1203 & ' can form a disulfide bridge?!!!'
1205 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1210 C Read preformed bridges.
1212 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1214 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1217 C Check if the residues involved in bridges are in the specified list of
1218 C bridging residues.
1221 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1222 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1223 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1224 & ' contains residues present in other pairs.'
1225 write (*,'(a,i3,a)') 'Disulfide pair',i,
1226 & ' contains residues present in other pairs.'
1228 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1234 if (ihpb(i).eq.iss(j)) goto 10
1236 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1239 if (jhpb(i).eq.iss(j)) goto 20
1241 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1247 ihpb(i)=ihpb(i)+nres
1248 jhpb(i)=jhpb(i)+nres
1254 c----------------------------------------------------------------------------
1255 subroutine read_x(kanal,*)
1256 implicit real*8 (a-h,o-z)
1257 include 'DIMENSIONS'
1258 include 'COMMON.GEO'
1259 include 'COMMON.VAR'
1260 include 'COMMON.CHAIN'
1261 include 'COMMON.IOUNITS'
1262 include 'COMMON.CONTROL'
1263 include 'COMMON.LOCAL'
1264 include 'COMMON.INTERACT'
1265 c Read coordinates from input
1267 read(kanal,'(8f10.5)',end=10,err=10)
1268 & ((c(l,k),l=1,3),k=1,nres),
1269 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1272 c(j,2*nres)=c(j,nres)
1274 call int_from_cart1(.false.)
1277 dc(j,i)=c(j,i+1)-c(j,i)
1278 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1282 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1284 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1285 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1293 c----------------------------------------------------------------------------
1294 subroutine read_threadbase
1295 implicit real*8 (a-h,o-z)
1296 include 'DIMENSIONS'
1297 include 'COMMON.IOUNITS'
1298 include 'COMMON.GEO'
1299 include 'COMMON.VAR'
1300 include 'COMMON.INTERACT'
1301 include 'COMMON.LOCAL'
1302 include 'COMMON.NAMES'
1303 include 'COMMON.CHAIN'
1304 include 'COMMON.FFIELD'
1305 include 'COMMON.SBRIDGE'
1306 include 'COMMON.HEADER'
1307 include 'COMMON.CONTROL'
1308 include 'COMMON.DBASE'
1309 include 'COMMON.THREAD'
1310 include 'COMMON.TIME1'
1311 C Read pattern database for threading.
1312 read (icbase,*) nseq
1314 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1315 & nres_base(2,i),nres_base(3,i)
1316 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1318 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1319 c & nres_base(2,i),nres_base(3,i)
1320 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1324 if (weidis.eq.0.0D0) weidis=0.1D0
1333 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1334 write (iout,'(a,i5)') 'nexcl: ',nexcl
1335 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1338 c------------------------------------------------------------------------------
1339 subroutine setup_var
1340 implicit real*8 (a-h,o-z)
1341 include 'DIMENSIONS'
1342 include 'COMMON.IOUNITS'
1343 include 'COMMON.GEO'
1344 include 'COMMON.VAR'
1345 include 'COMMON.INTERACT'
1346 include 'COMMON.LOCAL'
1347 include 'COMMON.NAMES'
1348 include 'COMMON.CHAIN'
1349 include 'COMMON.FFIELD'
1350 include 'COMMON.SBRIDGE'
1351 include 'COMMON.HEADER'
1352 include 'COMMON.CONTROL'
1353 include 'COMMON.DBASE'
1354 include 'COMMON.THREAD'
1355 include 'COMMON.TIME1'
1356 C Set up variable list.
1362 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1364 ialph(i,1)=nvar+nside
1368 if (indphi.gt.0) then
1370 else if (indback.gt.0) then
1375 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1378 c----------------------------------------------------------------------------
1379 subroutine gen_dist_constr
1380 C Generate CA distance constraints.
1381 implicit real*8 (a-h,o-z)
1382 include 'DIMENSIONS'
1383 include 'COMMON.IOUNITS'
1384 include 'COMMON.GEO'
1385 include 'COMMON.VAR'
1386 include 'COMMON.INTERACT'
1387 include 'COMMON.LOCAL'
1388 include 'COMMON.NAMES'
1389 include 'COMMON.CHAIN'
1390 include 'COMMON.FFIELD'
1391 include 'COMMON.SBRIDGE'
1392 include 'COMMON.HEADER'
1393 include 'COMMON.CONTROL'
1394 include 'COMMON.DBASE'
1395 include 'COMMON.THREAD'
1396 include 'COMMON.TIME1'
1397 dimension itype_pdb(maxres)
1398 common /pizda/ itype_pdb
1400 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1401 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1402 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1404 do i=nstart_sup,nstart_sup+nsup-1
1405 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1406 cd & ' seq_pdb', restyp(itype_pdb(i))
1407 do j=i+2,nstart_sup+nsup-1
1409 ihpb(nhpb)=i+nstart_seq-nstart_sup
1410 jhpb(nhpb)=j+nstart_seq-nstart_sup
1412 dhpb(nhpb)=dist(i,j)
1415 cd write (iout,'(a)') 'Distance constraints:'
1420 cd if (ii.gt.nres) then
1425 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1426 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1427 cd & dhpb(i),forcon(i)
1431 c----------------------------------------------------------------------------
1433 implicit real*8 (a-h,o-z)
1434 include 'DIMENSIONS'
1435 include 'COMMON.MAP'
1436 include 'COMMON.IOUNITS'
1437 character*3 angid(4) /'THE','PHI','ALP','OME'/
1438 character*80 mapcard,ucase
1440 read (inp,'(a)') mapcard
1441 mapcard=ucase(mapcard)
1442 if (index(mapcard,'PHI').gt.0) then
1444 else if (index(mapcard,'THE').gt.0) then
1446 else if (index(mapcard,'ALP').gt.0) then
1448 else if (index(mapcard,'OME').gt.0) then
1451 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1452 stop 'Error - illegal variable spec in MAP card.'
1454 call readi (mapcard,'RES1',res1(imap),0)
1455 call readi (mapcard,'RES2',res2(imap),0)
1456 if (res1(imap).eq.0) then
1457 res1(imap)=res2(imap)
1458 else if (res2(imap).eq.0) then
1459 res2(imap)=res1(imap)
1461 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1463 & 'Error - illegal definition of variable group in MAP.'
1464 stop 'Error - illegal definition of variable group in MAP.'
1466 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1467 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1468 call readi(mapcard,'NSTEP',nstep(imap),0)
1469 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1471 & 'Illegal boundary and/or step size specification in MAP.'
1472 stop 'Illegal boundary and/or step size specification in MAP.'
1477 c----------------------------------------------------------------------------
1479 implicit real*8 (a-h,o-z)
1480 include 'DIMENSIONS'
1481 include 'COMMON.IOUNITS'
1482 include 'COMMON.GEO'
1483 include 'COMMON.CSA'
1484 include 'COMMON.BANK'
1485 include 'COMMON.CONTROL'
1487 character*620 mcmcard
1488 call card_concat(mcmcard)
1490 call readi(mcmcard,'NCONF',nconf,50)
1491 call readi(mcmcard,'NADD',nadd,0)
1492 call readi(mcmcard,'JSTART',jstart,1)
1493 call readi(mcmcard,'JEND',jend,1)
1494 call readi(mcmcard,'NSTMAX',nstmax,500000)
1495 call readi(mcmcard,'N0',n0,1)
1496 call readi(mcmcard,'N1',n1,6)
1497 call readi(mcmcard,'N2',n2,4)
1498 call readi(mcmcard,'N3',n3,0)
1499 call readi(mcmcard,'N4',n4,0)
1500 call readi(mcmcard,'N5',n5,0)
1501 call readi(mcmcard,'N6',n6,10)
1502 call readi(mcmcard,'N7',n7,0)
1503 call readi(mcmcard,'N8',n8,0)
1504 call readi(mcmcard,'N9',n9,0)
1505 call readi(mcmcard,'N14',n14,0)
1506 call readi(mcmcard,'N15',n15,0)
1507 call readi(mcmcard,'N16',n16,0)
1508 call readi(mcmcard,'N17',n17,0)
1509 call readi(mcmcard,'N18',n18,0)
1511 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1513 call readi(mcmcard,'NDIFF',ndiff,2)
1514 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1515 call readi(mcmcard,'IS1',is1,1)
1516 call readi(mcmcard,'IS2',is2,8)
1517 call readi(mcmcard,'NRAN0',nran0,4)
1518 call readi(mcmcard,'NRAN1',nran1,2)
1519 call readi(mcmcard,'IRR',irr,1)
1520 call readi(mcmcard,'NSEED',nseed,20)
1521 call readi(mcmcard,'NTOTAL',ntotal,10000)
1522 call reada(mcmcard,'CUT1',cut1,2.0d0)
1523 call reada(mcmcard,'CUT2',cut2,5.0d0)
1524 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1525 call readi(mcmcard,'ICMAX',icmax,3)
1526 call readi(mcmcard,'IRESTART',irestart,0)
1527 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1530 call reada(mcmcard,'DELE',dele,20.0d0)
1531 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1532 call readi(mcmcard,'IREF',iref,0)
1533 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1534 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1535 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1536 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1537 write (iout,*) "NCONF_IN",nconf_in
1540 c----------------------------------------------------------------------------
1541 cfmc subroutine mcmfread
1542 cfmc implicit real*8 (a-h,o-z)
1543 cfmc include 'DIMENSIONS'
1544 cfmc include 'COMMON.MCMF'
1545 cfmc include 'COMMON.IOUNITS'
1546 cfmc include 'COMMON.GEO'
1547 cfmc character*80 ucase
1548 cfmc character*620 mcmcard
1549 cfmc call card_concat(mcmcard)
1551 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1552 cfmc write(iout,*)'MAXRANT=',maxrant
1553 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1554 cfmc write(iout,*)'MAXFAM=',maxfam
1555 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1556 cfmc write(iout,*)'NNET1=',nnet1
1557 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1558 cfmc write(iout,*)'NNET2=',nnet2
1559 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1560 cfmc write(iout,*)'NNET3=',nnet3
1561 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1562 cfmc write(iout,*)'ILASTT=',ilastt
1563 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1564 cfmc write(iout,*)'MAXSTR=',maxstr
1565 cfmc maxstr_f=maxstr/maxfam
1566 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1567 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1568 cfmc write(iout,*)'NMCMF=',nmcmf
1569 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1570 cfmc write(iout,*)'IFOCUS=',ifocus
1571 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1572 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1573 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1574 cfmc write(iout,*)'INTPRT=',intprt
1575 cfmc call readi(mcmcard,'IPRT',iprt,100)
1576 cfmc write(iout,*)'IPRT=',iprt
1577 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1578 cfmc write(iout,*)'IMAXTR=',imaxtr
1579 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1580 cfmc write(iout,*)'MAXEVEN=',maxeven
1581 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1582 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1583 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1584 cfmc write(iout,*)'INIMIN=',inimin
1585 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1586 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1587 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1588 cfmc write(iout,*)'NTHREAD=',nthread
1589 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1590 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1591 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1592 cfmc write(iout,*)'MAXPERT=',maxpert
1593 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1594 cfmc write(iout,*)'IRMSD=',irmsd
1595 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1596 cfmc write(iout,*)'DENEMIN=',denemin
1597 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1598 cfmc write(iout,*)'RCUT1S=',rcut1s
1599 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1600 cfmc write(iout,*)'RCUT1E=',rcut1e
1601 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1602 cfmc write(iout,*)'RCUT2S=',rcut2s
1603 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1604 cfmc write(iout,*)'RCUT2E=',rcut2e
1605 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1606 cfmc write(iout,*)'DPERT1=',d_pert1
1607 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1608 cfmc write(iout,*)'DPERT1A=',d_pert1a
1609 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1610 cfmc write(iout,*)'DPERT2=',d_pert2
1611 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1612 cfmc write(iout,*)'DPERT2A=',d_pert2a
1613 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1614 cfmc write(iout,*)'DPERT2B=',d_pert2b
1615 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1616 cfmc write(iout,*)'DPERT2C=',d_pert2c
1617 cfmc d_pert1=deg2rad*d_pert1
1618 cfmc d_pert1a=deg2rad*d_pert1a
1619 cfmc d_pert2=deg2rad*d_pert2
1620 cfmc d_pert2a=deg2rad*d_pert2a
1621 cfmc d_pert2b=deg2rad*d_pert2b
1622 cfmc d_pert2c=deg2rad*d_pert2c
1623 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1624 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1625 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1626 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1627 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1628 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1629 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1630 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1631 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1632 cfmc write(iout,*)'RCUTINI=',rcutini
1633 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1634 cfmc write(iout,*)'GRAT=',grat
1635 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1636 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1640 c----------------------------------------------------------------------------
1642 implicit real*8 (a-h,o-z)
1643 include 'DIMENSIONS'
1644 include 'COMMON.MCM'
1645 include 'COMMON.MCE'
1646 include 'COMMON.IOUNITS'
1648 character*320 mcmcard
1649 call card_concat(mcmcard)
1650 call readi(mcmcard,'MAXACC',maxacc,100)
1651 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1652 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1653 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1654 call readi(mcmcard,'MAXREPM',maxrepm,200)
1655 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1656 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1657 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1658 call reada(mcmcard,'E_UP',e_up,5.0D0)
1659 call reada(mcmcard,'DELTE',delte,0.1D0)
1660 call readi(mcmcard,'NSWEEP',nsweep,5)
1661 call readi(mcmcard,'NSTEPH',nsteph,0)
1662 call readi(mcmcard,'NSTEPC',nstepc,0)
1663 call reada(mcmcard,'TMIN',tmin,298.0D0)
1664 call reada(mcmcard,'TMAX',tmax,298.0D0)
1665 call readi(mcmcard,'NWINDOW',nwindow,0)
1666 call readi(mcmcard,'PRINT_MC',print_mc,0)
1667 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1668 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1669 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1670 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1671 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1672 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1673 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1674 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1675 if (nwindow.gt.0) then
1676 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1678 winlen(i)=winend(i)-winstart(i)+1
1681 if (tmax.lt.tmin) tmax=tmin
1682 if (tmax.eq.tmin) then
1686 if (nstepc.gt.0 .and. nsteph.gt.0) then
1687 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1688 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1690 C Probabilities of different move types
1691 sumpro_type(0)=0.0D0
1692 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1693 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1694 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1695 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1696 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1697 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1698 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1700 print *,'i',i,' sumprotype',sumpro_type(i)
1701 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1702 print *,'i',i,' sumprotype',sumpro_type(i)
1706 c----------------------------------------------------------------------------
1707 subroutine read_minim
1708 implicit real*8 (a-h,o-z)
1709 include 'DIMENSIONS'
1710 include 'COMMON.MINIM'
1711 include 'COMMON.IOUNITS'
1713 character*320 minimcard
1714 call card_concat(minimcard)
1715 call readi(minimcard,'MAXMIN',maxmin,2000)
1716 call readi(minimcard,'MAXFUN',maxfun,5000)
1717 call readi(minimcard,'MINMIN',minmin,maxmin)
1718 call readi(minimcard,'MINFUN',minfun,maxmin)
1719 call reada(minimcard,'TOLF',tolf,1.0D-2)
1720 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1721 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1722 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1723 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1724 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1725 & 'Options in energy minimization:'
1726 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1727 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1728 & 'MinMin:',MinMin,' MinFun:',MinFun,
1729 & ' TolF:',TolF,' RTolF:',RTolF
1732 c----------------------------------------------------------------------------
1733 subroutine read_angles(kanal,*)
1734 implicit real*8 (a-h,o-z)
1735 include 'DIMENSIONS'
1736 include 'COMMON.GEO'
1737 include 'COMMON.VAR'
1738 include 'COMMON.CHAIN'
1739 include 'COMMON.IOUNITS'
1740 include 'COMMON.CONTROL'
1741 c Read angles from input
1743 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1744 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1745 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1746 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1749 c 9/7/01 avoid 180 deg valence angle
1750 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1752 theta(i)=deg2rad*theta(i)
1753 phi(i)=deg2rad*phi(i)
1754 alph(i)=deg2rad*alph(i)
1755 omeg(i)=deg2rad*omeg(i)
1760 c----------------------------------------------------------------------------
1761 subroutine reada(rekord,lancuch,wartosc,default)
1763 character*(*) rekord,lancuch
1764 double precision wartosc,default
1767 iread=index(rekord,lancuch)
1768 if (iread.eq.0) then
1772 iread=iread+ilen(lancuch)+1
1773 read (rekord(iread:),*,err=10,end=10) wartosc
1778 c----------------------------------------------------------------------------
1779 subroutine readi(rekord,lancuch,wartosc,default)
1781 character*(*) rekord,lancuch
1782 integer wartosc,default
1785 iread=index(rekord,lancuch)
1786 if (iread.eq.0) then
1790 iread=iread+ilen(lancuch)+1
1791 read (rekord(iread:),*,err=10,end=10) wartosc
1796 c----------------------------------------------------------------------------
1797 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1800 integer tablica(dim),default
1801 character*(*) rekord,lancuch
1808 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1809 if (iread.eq.0) return
1810 iread=iread+ilen(lancuch)+1
1811 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1814 c----------------------------------------------------------------------------
1815 subroutine multreada(rekord,lancuch,tablica,dim,default)
1818 double precision tablica(dim),default
1819 character*(*) rekord,lancuch
1826 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1827 if (iread.eq.0) return
1828 iread=iread+ilen(lancuch)+1
1829 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1832 c----------------------------------------------------------------------------
1833 subroutine openunits
1834 implicit real*8 (a-h,o-z)
1835 include 'DIMENSIONS'
1838 character*16 form,nodename
1841 include 'COMMON.SETUP'
1842 include 'COMMON.IOUNITS'
1844 include 'COMMON.CONTROL'
1845 integer lenpre,lenpot,ilen,lentmp
1847 character*3 out1file_text,ucase
1850 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1851 call getenv_loc("PREFIX",prefix)
1853 call getenv_loc("POT",pot)
1854 call getenv_loc("DIRTMP",tmpdir)
1855 call getenv_loc("CURDIR",curdir)
1856 call getenv_loc("OUT1FILE",out1file_text)
1857 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1858 out1file_text=ucase(out1file_text)
1859 if (out1file_text(1:1).eq."Y") then
1862 out1file=fg_rank.gt.0
1867 if (lentmp.gt.0) then
1868 write (*,'(80(1h!))')
1869 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1870 write (*,'(80(1h!))')
1871 write (*,*)"All output files will be on node /tmp directory."
1873 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1874 if (me.eq.king) then
1875 write (*,*) "The master node is ",nodename
1876 else if (fg_rank.eq.0) then
1877 write (*,*) "I am the CG slave node ",nodename
1879 write (*,*) "I am the FG slave node ",nodename
1882 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1883 lenpre = lentmp+lenpre+1
1885 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1886 C Get the names and open the input files
1887 #if defined(WINIFL) || defined(WINPGI)
1888 open(1,file=pref_orig(:ilen(pref_orig))//
1889 & '.inp',status='old',readonly,shared)
1890 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1891 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1892 C Get parameter filenames and open the parameter files.
1893 call getenv_loc('BONDPAR',bondname)
1894 open (ibond,file=bondname,status='old',readonly,shared)
1895 call getenv_loc('THETPAR',thetname)
1896 open (ithep,file=thetname,status='old',readonly,shared)
1897 call getenv_loc('ROTPAR',rotname)
1898 open (irotam,file=rotname,status='old',readonly,shared)
1899 call getenv_loc('TORPAR',torname)
1900 open (itorp,file=torname,status='old',readonly,shared)
1901 call getenv_loc('TORDPAR',tordname)
1902 open (itordp,file=tordname,status='old',readonly,shared)
1903 call getenv_loc('FOURIER',fouriername)
1904 open (ifourier,file=fouriername,status='old',readonly,shared)
1905 call getenv_loc('ELEPAR',elename)
1906 open (ielep,file=elename,status='old',readonly,shared)
1907 call getenv_loc('SIDEPAR',sidename)
1908 open (isidep,file=sidename,status='old',readonly,shared)
1909 #elif (defined CRAY) || (defined AIX)
1910 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1912 c print *,"Processor",myrank," opened file 1"
1913 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1914 c print *,"Processor",myrank," opened file 9"
1915 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1916 C Get parameter filenames and open the parameter files.
1917 call getenv_loc('BONDPAR',bondname)
1918 open (ibond,file=bondname,status='old',action='read')
1919 c print *,"Processor",myrank," opened file IBOND"
1920 call getenv_loc('THETPAR',thetname)
1921 open (ithep,file=thetname,status='old',action='read')
1922 c print *,"Processor",myrank," opened file ITHEP"
1923 call getenv_loc('ROTPAR',rotname)
1924 open (irotam,file=rotname,status='old',action='read')
1925 c print *,"Processor",myrank," opened file IROTAM"
1926 call getenv_loc('TORPAR',torname)
1927 open (itorp,file=torname,status='old',action='read')
1928 c print *,"Processor",myrank," opened file ITORP"
1929 call getenv_loc('TORDPAR',tordname)
1930 open (itordp,file=tordname,status='old',action='read')
1931 c print *,"Processor",myrank," opened file ITORDP"
1932 call getenv_loc('SCCORPAR',sccorname)
1933 open (isccor,file=sccorname,status='old',action='read')
1934 c print *,"Processor",myrank," opened file ISCCOR"
1935 call getenv_loc('FOURIER',fouriername)
1936 open (ifourier,file=fouriername,status='old',action='read')
1937 c print *,"Processor",myrank," opened file IFOURIER"
1938 call getenv_loc('ELEPAR',elename)
1939 open (ielep,file=elename,status='old',action='read')
1940 c print *,"Processor",myrank," opened file IELEP"
1941 call getenv_loc('SIDEPAR',sidename)
1942 open (isidep,file=sidename,status='old',action='read')
1943 c print *,"Processor",myrank," opened file ISIDEP"
1944 c print *,"Processor",myrank," opened parameter files"
1946 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1947 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1948 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1949 C Get parameter filenames and open the parameter files.
1950 call getenv_loc('BONDPAR',bondname)
1951 open (ibond,file=bondname,status='old')
1952 call getenv_loc('THETPAR',thetname)
1953 open (ithep,file=thetname,status='old')
1954 call getenv_loc('ROTPAR',rotname)
1955 open (irotam,file=rotname,status='old')
1956 call getenv_loc('TORPAR',torname)
1957 open (itorp,file=torname,status='old')
1958 call getenv_loc('TORDPAR',tordname)
1959 open (itordp,file=tordname,status='old')
1960 call getenv_loc('SCCORPAR',sccorname)
1961 open (isccor,file=sccorname,status='old')
1962 call getenv_loc('FOURIER',fouriername)
1963 open (ifourier,file=fouriername,status='old')
1964 call getenv_loc('ELEPAR',elename)
1965 open (ielep,file=elename,status='old')
1966 call getenv_loc('SIDEPAR',sidename)
1967 open (isidep,file=sidename,status='old')
1969 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1971 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1972 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1973 C Get parameter filenames and open the parameter files.
1974 call getenv_loc('BONDPAR',bondname)
1975 open (ibond,file=bondname,status='old',readonly)
1976 call getenv_loc('THETPAR',thetname)
1977 open (ithep,file=thetname,status='old',readonly)
1978 call getenv_loc('ROTPAR',rotname)
1979 open (irotam,file=rotname,status='old',readonly)
1980 call getenv_loc('TORPAR',torname)
1981 open (itorp,file=torname,status='old',readonly)
1982 call getenv_loc('TORDPAR',tordname)
1983 open (itordp,file=tordname,status='old',readonly)
1984 call getenv_loc('SCCORPAR',sccorname)
1985 open (isccor,file=sccorname,status='old',readonly)
1987 call getenv_loc('THETPARPDB',thetname_pdb)
1988 print *,"thetname_pdb ",thetname_pdb
1989 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1990 print *,ithep_pdb," opened"
1992 call getenv_loc('FOURIER',fouriername)
1993 open (ifourier,file=fouriername,status='old',readonly)
1994 call getenv_loc('ELEPAR',elename)
1995 open (ielep,file=elename,status='old',readonly)
1996 call getenv_loc('SIDEPAR',sidename)
1997 open (isidep,file=sidename,status='old',readonly)
1998 call getenv_loc('LIPTRANPAR',liptranname)
1999 open (iliptranpar,file=liptranname,status='old',action='read')
2001 call getenv_loc('ROTPARPDB',rotname_pdb)
2002 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2007 C 8/9/01 In the newest version SCp interaction constants are read from a file
2008 C Use -DOLDSCP to use hard-coded constants instead.
2010 call getenv_loc('SCPPAR',scpname)
2011 #if defined(WINIFL) || defined(WINPGI)
2012 open (iscpp,file=scpname,status='old',readonly,shared)
2013 #elif (defined CRAY) || (defined AIX)
2014 open (iscpp,file=scpname,status='old',action='read')
2016 open (iscpp,file=scpname,status='old')
2018 open (iscpp,file=scpname,status='old',readonly)
2021 call getenv_loc('PATTERN',patname)
2022 #if defined(WINIFL) || defined(WINPGI)
2023 open (icbase,file=patname,status='old',readonly,shared)
2024 #elif (defined CRAY) || (defined AIX)
2025 open (icbase,file=patname,status='old',action='read')
2027 open (icbase,file=patname,status='old')
2029 open (icbase,file=patname,status='old',readonly)
2032 C Open output file only for CG processes
2033 c print *,"Processor",myrank," fg_rank",fg_rank
2034 if (fg_rank.eq.0) then
2036 if (nodes.eq.1) then
2039 npos = dlog10(dfloat(nodes-1))+1
2041 if (npos.lt.3) npos=3
2042 write (liczba,'(i1)') npos
2043 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2045 write (liczba,form) me
2046 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2047 & liczba(:ilen(liczba))
2048 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2050 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2052 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2053 & liczba(:ilen(liczba))//'.mol2'
2054 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2055 & liczba(:ilen(liczba))//'.stat'
2057 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2058 & //liczba(:ilen(liczba))//'.stat')
2059 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2062 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2063 & liczba(:ilen(liczba))//'.const'
2068 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2069 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2070 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2071 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2072 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2074 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2076 rest2name=prefix(:ilen(prefix))//'.rst'
2078 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2081 #if defined(AIX) || defined(PGI)
2082 if (me.eq.king .or. .not. out1file)
2083 & open(iout,file=outname,status='unknown')
2085 if (fg_rank.gt.0) then
2086 write (liczba,'(i3.3)') myrank/nfgtasks
2087 write (ll,'(bz,i3.3)') fg_rank
2088 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2093 open(igeom,file=intname,status='unknown',position='append')
2094 open(ipdb,file=pdbname,status='unknown')
2095 open(imol2,file=mol2name,status='unknown')
2096 open(istat,file=statname,status='unknown',position='append')
2098 c1out open(iout,file=outname,status='unknown')
2101 if (me.eq.king .or. .not.out1file)
2102 & open(iout,file=outname,status='unknown')
2104 if (fg_rank.gt.0) then
2105 write (liczba,'(i3.3)') myrank/nfgtasks
2106 write (ll,'(bz,i3.3)') fg_rank
2107 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2112 open(igeom,file=intname,status='unknown',access='append')
2113 open(ipdb,file=pdbname,status='unknown')
2114 open(imol2,file=mol2name,status='unknown')
2115 open(istat,file=statname,status='unknown',access='append')
2117 c1out open(iout,file=outname,status='unknown')
2120 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2121 csa_seed=prefix(:lenpre)//'.CSA.seed'
2122 csa_history=prefix(:lenpre)//'.CSA.history'
2123 csa_bank=prefix(:lenpre)//'.CSA.bank'
2124 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2125 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2126 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2127 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2128 csa_int=prefix(:lenpre)//'.int'
2129 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2130 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2131 csa_in=prefix(:lenpre)//'.CSA.in'
2132 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2135 write (iout,'(80(1h-))')
2136 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2137 write (iout,'(80(1h-))')
2138 write (iout,*) "Input file : ",
2139 & pref_orig(:ilen(pref_orig))//'.inp'
2140 write (iout,*) "Output file : ",
2141 & outname(:ilen(outname))
2143 write (iout,*) "Sidechain potential file : ",
2144 & sidename(:ilen(sidename))
2146 write (iout,*) "SCp potential file : ",
2147 & scpname(:ilen(scpname))
2149 write (iout,*) "Electrostatic potential file : ",
2150 & elename(:ilen(elename))
2151 write (iout,*) "Cumulant coefficient file : ",
2152 & fouriername(:ilen(fouriername))
2153 write (iout,*) "Torsional parameter file : ",
2154 & torname(:ilen(torname))
2155 write (iout,*) "Double torsional parameter file : ",
2156 & tordname(:ilen(tordname))
2157 write (iout,*) "SCCOR parameter file : ",
2158 & sccorname(:ilen(sccorname))
2159 write (iout,*) "Bond & inertia constant file : ",
2160 & bondname(:ilen(bondname))
2161 write (iout,*) "Bending parameter file : ",
2162 & thetname(:ilen(thetname))
2163 write (iout,*) "Rotamer parameter file : ",
2164 & rotname(:ilen(rotname))
2165 write (iout,*) "Thetpdb parameter file : ",
2166 & thetname_pdb(:ilen(thetname_pdb))
2167 write (iout,*) "Threading database : ",
2168 & patname(:ilen(patname))
2170 &write (iout,*)" DIRTMP : ",
2172 write (iout,'(80(1h-))')
2176 c----------------------------------------------------------------------------
2177 subroutine card_concat(card)
2178 implicit real*8 (a-h,o-z)
2179 include 'DIMENSIONS'
2180 include 'COMMON.IOUNITS'
2182 character*80 karta,ucase
2184 read (inp,'(a)') karta
2187 do while (karta(80:80).eq.'&')
2188 card=card(:ilen(card)+1)//karta(:79)
2189 read (inp,'(a)') karta
2192 card=card(:ilen(card)+1)//karta
2195 c----------------------------------------------------------------------------------
2197 implicit real*8 (a-h,o-z)
2198 include 'DIMENSIONS'
2199 include 'COMMON.CHAIN'
2200 include 'COMMON.IOUNITS'
2202 open(irest2,file=rest2name,status='unknown')
2203 read(irest2,*) totT,EK,potE,totE,t_bath
2205 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2208 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2211 read (irest2,*) iset
2216 c---------------------------------------------------------------------------------
2217 subroutine read_fragments
2218 implicit real*8 (a-h,o-z)
2219 include 'DIMENSIONS'
2223 include 'COMMON.SETUP'
2224 include 'COMMON.CHAIN'
2225 include 'COMMON.IOUNITS'
2227 include 'COMMON.CONTROL'
2228 read(inp,*) nset,nfrag,npair,nfrag_back
2229 if(me.eq.king.or..not.out1file)
2230 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2231 & " nfrag_back",nfrag_back
2233 read(inp,*) mset(iset)
2235 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2237 if(me.eq.king.or..not.out1file)
2238 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2239 & ifrag(2,i,iset), qinfrag(i,iset)
2242 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2244 if(me.eq.king.or..not.out1file)
2245 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2246 & ipair(2,i,iset), qinpair(i,iset)
2249 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2250 & wfrag_back(3,i,iset),
2251 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2252 if(me.eq.king.or..not.out1file)
2253 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2254 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2259 C---------------------------------------------------------------------------
2260 subroutine read_afminp
2261 implicit real*8 (a-h,o-z)
2262 include 'DIMENSIONS'
2266 include 'COMMON.SETUP'
2267 include 'COMMON.CONTROL'
2268 include 'COMMON.CHAIN'
2269 include 'COMMON.IOUNITS'
2270 include 'COMMON.SBRIDGE'
2271 character*320 afmcard
2273 call card_concat(afmcard)
2274 call readi(afmcard,"BEG",afmbeg,0)
2275 call readi(afmcard,"END",afmend,0)
2276 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2277 print *,'FORCE=' ,forceAFMconst
2278 CCCC NOW PROPERTIES FOR AFM
2281 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2283 distafminit=dsqrt(distafminit)
2284 print *,'initdist',distafminit
2288 c-------------------------------------------------------------------------------
2289 subroutine read_dist_constr
2290 implicit real*8 (a-h,o-z)
2291 include 'DIMENSIONS'
2295 include 'COMMON.SETUP'
2296 include 'COMMON.CONTROL'
2297 include 'COMMON.CHAIN'
2298 include 'COMMON.IOUNITS'
2299 include 'COMMON.SBRIDGE'
2300 integer ifrag_(2,100),ipair_(2,100)
2301 double precision wfrag_(100),wpair_(100)
2302 character*500 controlcard
2303 c write (iout,*) "Calling read_dist_constr"
2304 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2306 call card_concat(controlcard)
2307 call readi(controlcard,"NFRAG",nfrag_,0)
2308 call readi(controlcard,"NPAIR",npair_,0)
2309 call readi(controlcard,"NDIST",ndist_,0)
2310 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2311 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2312 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2313 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2314 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2315 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2316 c write (iout,*) "IFRAG"
2318 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2320 c write (iout,*) "IPAIR"
2322 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2326 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2327 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2328 & ifrag_(2,i)=nstart_sup+nsup-1
2329 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2331 if (wfrag_(i).gt.0.0d0) then
2332 do j=ifrag_(1,i),ifrag_(2,i)-1
2333 do k=j+1,ifrag_(2,i)
2334 c write (iout,*) "j",j," k",k
2336 if (constr_dist.eq.1) then
2341 forcon(nhpb)=wfrag_(i)
2342 else if (constr_dist.eq.2) then
2343 if (ddjk.le.dist_cut) then
2348 forcon(nhpb)=wfrag_(i)
2355 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2358 if (.not.out1file .or. me.eq.king)
2359 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2360 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2362 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2363 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2370 if (wpair_(i).gt.0.0d0) then
2378 do j=ifrag_(1,ii),ifrag_(2,ii)
2379 do k=ifrag_(1,jj),ifrag_(2,jj)
2383 forcon(nhpb)=wpair_(i)
2384 dhpb(nhpb)=dist(j,k)
2386 if (.not.out1file .or. me.eq.king)
2387 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2388 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2390 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2391 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2398 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2399 if (forcon(nhpb+1).gt.0.0d0) then
2401 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2403 if (.not.out1file .or. me.eq.king)
2404 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2405 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2407 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2408 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2415 c-------------------------------------------------------------------------------
2417 subroutine flush(iu)
2422 subroutine flush(iu)
2427 c------------------------------------------------------------------------------
2428 subroutine copy_to_tmp(source)
2429 include "DIMENSIONS"
2430 include "COMMON.IOUNITS"
2431 character*(*) source
2432 character* 256 tmpfile
2436 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2437 inquire(file=tmpfile,exist=ex)
2439 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2440 & " to temporary directory..."
2441 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2442 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2446 c------------------------------------------------------------------------------
2447 subroutine move_from_tmp(source)
2448 include "DIMENSIONS"
2449 include "COMMON.IOUNITS"
2450 character*(*) source
2453 write (*,*) "Moving ",source(:ilen(source)),
2454 & " from temporary directory to working directory"
2455 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2456 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2459 c------------------------------------------------------------------------------
2460 subroutine random_init(seed)
2462 C Initialize random number generator
2464 implicit real*8 (a-h,o-z)
2465 include 'DIMENSIONS'
2468 logical OKRandom, prng_restart
2470 integer iseed_array(4)
2472 include 'COMMON.IOUNITS'
2473 include 'COMMON.TIME1'
2474 include 'COMMON.THREAD'
2475 include 'COMMON.SBRIDGE'
2476 include 'COMMON.CONTROL'
2477 include 'COMMON.MCM'
2478 include 'COMMON.MAP'
2479 include 'COMMON.HEADER'
2480 include 'COMMON.CSA'
2481 include 'COMMON.CHAIN'
2482 include 'COMMON.MUCA'
2484 include 'COMMON.FFIELD'
2485 include 'COMMON.SETUP'
2486 iseed=-dint(dabs(seed))
2487 if (iseed.eq.0) then
2488 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2489 & 'Random seed undefined. The program will stop.'
2490 write (*,'(/80(1h*)/20x,a/80(1h*))')
2491 & 'Random seed undefined. The program will stop.'
2493 call mpi_finalize(mpi_comm_world,ierr)
2495 stop 'Bad random seed.'
2498 if (fg_rank.eq.0) then
2502 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2503 OKRandom = prng_restart(me,iseed)
2506 tmp=65536.0d0**(4-i)
2507 iseed_array(i) = dint(seed/tmp)
2508 seed=seed-iseed_array(i)*tmp
2511 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2512 & (iseed_array(i),i=1,4)
2513 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2514 & (iseed_array(i),i=1,4)
2515 OKRandom = prng_restart(me,iseed_array)
2518 c r1 = prng_next(me)
2519 r1=ran_number(0.0D0,1.0D0)
2521 & write (iout,*) 'ran_num',r1
2522 if (r1.lt.0.0d0) OKRandom=.false.
2524 if (.not.OKRandom) then
2525 write (iout,*) 'PRNG IS NOT WORKING!!!'
2526 print *,'PRNG IS NOT WORKING!!!'
2529 call mpi_abort(mpi_comm_world,error_msg,ierr)
2532 write (iout,*) 'too many processors for parallel prng'
2533 write (*,*) 'too many processors for parallel prng'
2541 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)