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 selfguide=(index(controlcard,'SELFGUIDE'))
142 print *,'AFMlog',AFMlog,selfguide,"KUPA"
143 call readi(controlcard,'IPRINT',iprint,0)
144 call readi(controlcard,'MAXGEN',maxgen,10000)
145 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
146 call readi(controlcard,"KDIAG",kdiag,0)
147 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
148 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
149 & write (iout,*) "RESCALE_MODE",rescale_mode
150 split_ene=index(controlcard,'SPLIT_ENE').gt.0
151 if (index(controlcard,'REGULAR').gt.0.0D0) then
152 call reada(controlcard,'WEIDIS',weidis,0.1D0)
156 if (index(controlcard,'CHECKGRAD').gt.0) then
158 if (index(controlcard,'CART').gt.0) then
160 elseif (index(controlcard,'CARINT').gt.0) then
165 elseif (index(controlcard,'THREAD').gt.0) then
167 call readi(controlcard,'THREAD',nthread,0)
168 if (nthread.gt.0) then
169 call reada(controlcard,'WEIDIS',weidis,0.1D0)
172 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
173 stop 'Error termination in Read_Control.'
175 else if (index(controlcard,'MCMA').gt.0) then
177 else if (index(controlcard,'MCEE').gt.0) then
179 else if (index(controlcard,'MULTCONF').gt.0) then
181 else if (index(controlcard,'MAP').gt.0) then
183 call readi(controlcard,'MAP',nmap,0)
184 else if (index(controlcard,'CSA').gt.0) then
186 crc else if (index(controlcard,'ZSCORE').gt.0) then
188 crc ZSCORE is rm from UNRES, modecalc=9 is available
191 cfcm else if (index(controlcard,'MCMF').gt.0) then
193 else if (index(controlcard,'SOFTREG').gt.0) then
195 else if (index(controlcard,'CHECK_BOND').gt.0) then
197 else if (index(controlcard,'TEST').gt.0) then
199 else if (index(controlcard,'MD').gt.0) then
201 else if (index(controlcard,'RE ').gt.0) then
205 lmuca=index(controlcard,'MUCA').gt.0
206 call readi(controlcard,'MUCADYN',mucadyn,0)
207 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
208 if (lmuca .and. (me.eq.king .or. .not.out1file ))
210 write (iout,*) 'MUCADYN=',mucadyn
211 write (iout,*) 'MUCASMOOTH=',muca_smooth
214 iscode=index(controlcard,'ONE_LETTER')
215 indphi=index(controlcard,'PHI')
216 indback=index(controlcard,'BACK')
217 iranconf=index(controlcard,'RAND_CONF')
218 i2ndstr=index(controlcard,'USE_SEC_PRED')
219 gradout=index(controlcard,'GRADOUT').gt.0
220 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
221 C DISTCHAINMAX become obsolete for periodic boundry condition
222 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
223 C Reading the dimensions of box in x,y,z coordinates
224 call reada(controlcard,'BOXX',boxxsize,100.0d0)
225 call reada(controlcard,'BOXY',boxysize,100.0d0)
226 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
227 c Cutoff range for interactions
228 call reada(controlcard,"R_CUT",r_cut,15.0d0)
229 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
230 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
231 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
232 if (lipthick.gt.0.0d0) then
233 bordliptop=(boxzsize+lipthick)/2.0
234 bordlipbot=bordliptop-lipthick
236 if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0))
237 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
238 buflipbot=bordlipbot+lipbufthick
239 bufliptop=bordliptop-lipbufthick
240 if ((lipbufthick*2.0d0).gt.lipthick)
241 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
243 write(iout,*) "bordliptop=",bordliptop
244 write(iout,*) "bordlipbot=",bordlipbot
245 write(iout,*) "bufliptop=",bufliptop
246 write(iout,*) "buflipbot=",buflipbot
249 if (me.eq.king .or. .not.out1file )
250 & write (iout,*) "DISTCHAINMAX",distchainmax
252 if(me.eq.king.or..not.out1file)
253 & write (iout,'(2a)') diagmeth(kdiag),
254 & ' routine used to diagonalize matrices.'
257 c--------------------------------------------------------------------------
258 subroutine read_REMDpar
262 implicit real*8 (a-h,o-z)
264 include 'COMMON.IOUNITS'
265 include 'COMMON.TIME1'
268 include 'COMMON.LANGEVIN'
270 include 'COMMON.LANGEVIN.lang0'
272 include 'COMMON.INTERACT'
273 include 'COMMON.NAMES'
275 include 'COMMON.REMD'
276 include 'COMMON.CONTROL'
277 include 'COMMON.SETUP'
279 character*320 controlcard
280 character*3200 controlcard1
281 integer iremd_m_total
283 if(me.eq.king.or..not.out1file)
284 & write (iout,*) "REMD setup"
286 call card_concat(controlcard)
287 call readi(controlcard,"NREP",nrep,3)
288 call readi(controlcard,"NSTEX",nstex,1000)
289 call reada(controlcard,"RETMIN",retmin,10.0d0)
290 call reada(controlcard,"RETMAX",retmax,1000.0d0)
291 mremdsync=(index(controlcard,'SYNC').gt.0)
292 call readi(controlcard,"NSYN",i_sync_step,100)
293 restart1file=(index(controlcard,'REST1FILE').gt.0)
294 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
295 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
296 if(max_cache_traj_use.gt.max_cache_traj)
297 & max_cache_traj_use=max_cache_traj
298 if(me.eq.king.or..not.out1file) then
299 cd if (traj1file) then
300 crc caching is in testing - NTWX is not ignored
301 cd write (iout,*) "NTWX value is ignored"
302 cd write (iout,*) " trajectory is stored to one file by master"
303 cd write (iout,*) " before exchange at NSTEX intervals"
305 write (iout,*) "NREP= ",nrep
306 write (iout,*) "NSTEX= ",nstex
307 write (iout,*) "SYNC= ",mremdsync
308 write (iout,*) "NSYN= ",i_sync_step
309 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
312 if (index(controlcard,'TLIST').gt.0) then
314 call card_concat(controlcard1)
315 read(controlcard1,*) (remd_t(i),i=1,nrep)
316 if(me.eq.king.or..not.out1file)
317 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
320 if (index(controlcard,'MLIST').gt.0) then
322 call card_concat(controlcard1)
323 read(controlcard1,*) (remd_m(i),i=1,nrep)
324 if(me.eq.king.or..not.out1file) then
325 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
328 iremd_m_total=iremd_m_total+remd_m(i)
330 write (iout,*) 'Total number of replicas ',iremd_m_total
333 if(me.eq.king.or..not.out1file)
334 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
337 c--------------------------------------------------------------------------
338 subroutine read_MDpar
342 implicit real*8 (a-h,o-z)
344 include 'COMMON.IOUNITS'
345 include 'COMMON.TIME1'
348 include 'COMMON.LANGEVIN'
350 include 'COMMON.LANGEVIN.lang0'
352 include 'COMMON.INTERACT'
353 include 'COMMON.NAMES'
355 include 'COMMON.SETUP'
356 include 'COMMON.CONTROL'
357 include 'COMMON.SPLITELE'
359 character*320 controlcard
361 call card_concat(controlcard)
362 call readi(controlcard,"NSTEP",n_timestep,1000000)
363 call readi(controlcard,"NTWE",ntwe,100)
364 call readi(controlcard,"NTWX",ntwx,1000)
365 call reada(controlcard,"DT",d_time,1.0d-1)
366 call reada(controlcard,"DVMAX",dvmax,2.0d1)
367 call reada(controlcard,"DAMAX",damax,1.0d1)
368 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
369 call readi(controlcard,"LANG",lang,0)
370 RESPA = index(controlcard,"RESPA") .gt. 0
371 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
372 ntime_split0=ntime_split
373 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
374 ntime_split0=ntime_split
375 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
376 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
377 rest = index(controlcard,"REST").gt.0
378 tbf = index(controlcard,"TBF").gt.0
379 usampl = index(controlcard,"USAMPL").gt.0
381 mdpdb = index(controlcard,"MDPDB").gt.0
382 call reada(controlcard,"T_BATH",t_bath,300.0d0)
383 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
384 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
385 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
386 if (count_reset_moment.eq.0) count_reset_moment=1000000000
387 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
388 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
389 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
390 if (count_reset_vel.eq.0) count_reset_vel=1000000000
391 large = index(controlcard,"LARGE").gt.0
392 print_compon = index(controlcard,"PRINT_COMPON").gt.0
393 rattle = index(controlcard,"RATTLE").gt.0
394 c if performing umbrella sampling, fragments constrained are read from the fragment file
400 if(me.eq.king.or..not.out1file) then
402 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
404 write (iout,'(a)') "The units are:"
405 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
406 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
407 & " acceleration: angstrom/(48.9 fs)**2"
408 write (iout,'(a)') "energy: kcal/mol, temperature: K"
410 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
411 write (iout,'(a60,f10.5,a)')
412 & "Initial time step of numerical integration:",d_time,
414 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
416 write (iout,'(2a,i4,a)')
417 & "A-MTS algorithm used; initial time step for fast-varying",
418 & " short-range forces split into",ntime_split," steps."
419 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
420 & r_cut," lambda",rlamb
422 write (iout,'(2a,f10.5)')
423 & "Maximum acceleration threshold to reduce the time step",
424 & "/increase split number:",damax
425 write (iout,'(2a,f10.5)')
426 & "Maximum predicted energy drift to reduce the timestep",
427 & "/increase split number:",edriftmax
428 write (iout,'(a60,f10.5)')
429 & "Maximum velocity threshold to reduce velocities:",dvmax
430 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
431 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
432 if (rattle) write (iout,'(a60)')
433 & "Rattle algorithm used to constrain the virtual bonds"
437 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
438 call reada(controlcard,"RWAT",rwat,1.4d0)
439 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
440 surfarea=index(controlcard,"SURFAREA").gt.0
441 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
442 if(me.eq.king.or..not.out1file)then
443 write (iout,'(/a,$)') "Langevin dynamics calculation"
446 & " with direct integration of Langevin equations"
447 else if (lang.eq.2) then
448 write (iout,'(a/)') " with TINKER stochasic MD integrator"
449 else if (lang.eq.3) then
450 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
451 else if (lang.eq.4) then
452 write (iout,'(a/)') " in overdamped mode"
454 write (iout,'(//a,i5)')
455 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
458 write (iout,'(a60,f10.5)') "Temperature:",t_bath
459 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
460 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
461 write (iout,'(a60,f10.5)')
462 & "Scaling factor of the friction forces:",scal_fric
463 if (surfarea) write (iout,'(2a,i10,a)')
464 & "Friction coefficients will be scaled by solvent-accessible",
465 & " surface area every",reset_fricmat," steps."
467 c Calculate friction coefficients and bounds of stochastic forces
468 eta=6*pi*cPoise*etawat
469 if(me.eq.king.or..not.out1file)
470 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
472 gamp=scal_fric*(pstok+rwat)*eta
473 stdfp=dsqrt(2*Rb*t_bath/d_time)
475 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
476 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
478 if(me.eq.king.or..not.out1file)then
479 write (iout,'(/2a/)')
480 & "Radii of site types and friction coefficients and std's of",
481 & " stochastic forces of fully exposed sites"
482 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
484 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
485 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
489 if(me.eq.king.or..not.out1file)then
490 write (iout,'(a)') "Berendsen bath calculation"
491 write (iout,'(a60,f10.5)') "Temperature:",t_bath
492 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
494 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
495 & count_reset_moment," steps"
497 & write (iout,'(a,i10,a)')
498 & "Velocities will be reset at random every",count_reset_vel,
502 if(me.eq.king.or..not.out1file)
503 & write (iout,'(a31)') "Microcanonical mode calculation"
505 if(me.eq.king.or..not.out1file)then
506 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
508 write(iout,*) "MD running with constraints."
509 write(iout,*) "Equilibration time ", eq_time, " mtus."
510 write(iout,*) "Constraining ", nfrag," fragments."
511 write(iout,*) "Length of each fragment, weight and q0:"
513 write (iout,*) "Set of restraints #",iset
515 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
516 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
518 write(iout,*) "constraints between ", npair, "fragments."
519 write(iout,*) "constraint pairs, weights and q0:"
521 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
522 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
524 write(iout,*) "angle constraints within ", nfrag_back,
525 & "backbone fragments."
526 write(iout,*) "fragment, weights:"
528 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
529 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
530 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
533 iset=mod(kolor,nset)+1
536 if(me.eq.king.or..not.out1file)
537 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
540 c------------------------------------------------------------------------------
543 C Read molecular data.
545 implicit real*8 (a-h,o-z)
551 include 'COMMON.IOUNITS'
554 include 'COMMON.INTERACT'
555 include 'COMMON.LOCAL'
556 include 'COMMON.NAMES'
557 include 'COMMON.CHAIN'
558 include 'COMMON.FFIELD'
559 include 'COMMON.SBRIDGE'
560 include 'COMMON.HEADER'
561 include 'COMMON.CONTROL'
562 include 'COMMON.DBASE'
563 include 'COMMON.THREAD'
564 include 'COMMON.CONTACTS'
565 include 'COMMON.TORCNSTR'
566 include 'COMMON.TIME1'
567 include 'COMMON.BOUNDS'
569 include 'COMMON.SETUP'
570 character*4 sequence(maxres)
572 double precision x(maxvar)
573 character*256 pdbfile
574 character*320 weightcard
575 character*80 weightcard_t,ucase
576 dimension itype_pdb(maxres)
577 common /pizda/ itype_pdb
578 logical seq_comp,fail
579 double precision energia(0:n_ene)
585 C Read weights of the subsequent energy terms.
586 call card_concat(weightcard)
587 call reada(weightcard,'WLONG',wlong,1.0D0)
588 call reada(weightcard,'WSC',wsc,wlong)
589 call reada(weightcard,'WSCP',wscp,wlong)
590 call reada(weightcard,'WELEC',welec,1.0D0)
591 call reada(weightcard,'WVDWPP',wvdwpp,welec)
592 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
593 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
594 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
595 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
596 call reada(weightcard,'WTURN3',wturn3,1.0D0)
597 call reada(weightcard,'WTURN4',wturn4,1.0D0)
598 call reada(weightcard,'WTURN6',wturn6,1.0D0)
599 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
600 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
601 call reada(weightcard,'WBOND',wbond,1.0D0)
602 call reada(weightcard,'WTOR',wtor,1.0D0)
603 call reada(weightcard,'WTORD',wtor_d,1.0D0)
604 call reada(weightcard,'WANG',wang,1.0D0)
605 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
606 call reada(weightcard,'SCAL14',scal14,0.4D0)
607 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
608 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
609 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
610 call reada(weightcard,'TEMP0',temp0,300.0d0)
611 call reada(weightcard,'WLT',wliptran,0.0D0)
612 if (index(weightcard,'SOFT').gt.0) ipot=6
613 C 12/1/95 Added weight for the multi-body term WCORR
614 call reada(weightcard,'WCORRH',wcorr,1.0D0)
615 if (wcorr4.gt.0.0d0) wcorr=wcorr4
635 if(me.eq.king.or..not.out1file)
636 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
637 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
639 10 format (/'Energy-term weights (unscaled):'//
640 & 'WSCC= ',f10.6,' (SC-SC)'/
641 & 'WSCP= ',f10.6,' (SC-p)'/
642 & 'WELEC= ',f10.6,' (p-p electr)'/
643 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
644 & 'WBOND= ',f10.6,' (stretching)'/
645 & 'WANG= ',f10.6,' (bending)'/
646 & 'WSCLOC= ',f10.6,' (SC local)'/
647 & 'WTOR= ',f10.6,' (torsional)'/
648 & 'WTORD= ',f10.6,' (double torsional)'/
649 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
650 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
651 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
652 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
653 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
654 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
655 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
656 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
657 & 'WTURN6= ',f10.6,' (turns, 6th order)')
658 if(me.eq.king.or..not.out1file)then
659 if (wcorr4.gt.0.0d0) then
660 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
661 & 'between contact pairs of peptide groups'
662 write (iout,'(2(a,f5.3/))')
663 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
664 & 'Range of quenching the correlation terms:',2*delt_corr
665 else if (wcorr.gt.0.0d0) then
666 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
667 & 'between contact pairs of peptide groups'
669 write (iout,'(a,f8.3)')
670 & 'Scaling factor of 1,4 SC-p interactions:',scal14
671 write (iout,'(a,f8.3)')
672 & 'General scaling factor of SC-p interactions:',scalscp
674 r0_corr=cutoff_corr-delt_corr
676 aad(i,1)=scalscp*aad(i,1)
677 aad(i,2)=scalscp*aad(i,2)
678 bad(i,1)=scalscp*bad(i,1)
679 bad(i,2)=scalscp*bad(i,2)
681 call rescale_weights(t_bath)
682 if(me.eq.king.or..not.out1file)
683 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
684 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
686 22 format (/'Energy-term weights (scaled):'//
687 & 'WSCC= ',f10.6,' (SC-SC)'/
688 & 'WSCP= ',f10.6,' (SC-p)'/
689 & 'WELEC= ',f10.6,' (p-p electr)'/
690 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
691 & 'WBOND= ',f10.6,' (stretching)'/
692 & 'WANG= ',f10.6,' (bending)'/
693 & 'WSCLOC= ',f10.6,' (SC local)'/
694 & 'WTOR= ',f10.6,' (torsional)'/
695 & 'WTORD= ',f10.6,' (double torsional)'/
696 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
697 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
698 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
699 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
700 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
701 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
702 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
703 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
704 & 'WTURN6= ',f10.6,' (turns, 6th order)')
705 if(me.eq.king.or..not.out1file)
706 & write (iout,*) "Reference temperature for weights calculation:",
708 call reada(weightcard,"D0CM",d0cm,3.78d0)
709 call reada(weightcard,"AKCM",akcm,15.1d0)
710 call reada(weightcard,"AKTH",akth,11.0d0)
711 call reada(weightcard,"AKCT",akct,12.0d0)
712 call reada(weightcard,"V1SS",v1ss,-1.08d0)
713 call reada(weightcard,"V2SS",v2ss,7.61d0)
714 call reada(weightcard,"V3SS",v3ss,13.7d0)
715 call reada(weightcard,"EBR",ebr,-5.50D0)
716 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
718 dyn_ss_mask(i)=.false.
722 dyn_ssbond_ij(i,j)=1.0d300
725 call reada(weightcard,"HT",Ht,0.0D0)
727 ss_depth=ebr/wsc-0.25*eps(1,1)
728 Ht=Ht/wsc-0.25*eps(1,1)
729 akcm=akcm*wstrain/wsc
730 akth=akth*wstrain/wsc
731 akct=akct*wstrain/wsc
732 v1ss=v1ss*wstrain/wsc
733 v2ss=v2ss*wstrain/wsc
734 v3ss=v3ss*wstrain/wsc
736 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
739 if(me.eq.king.or..not.out1file) then
740 write (iout,*) "Parameters of the SS-bond potential:"
741 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
743 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
744 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
745 write (iout,*)" HT",Ht
746 print *,'indpdb=',indpdb,' pdbref=',pdbref
748 if (indpdb.gt.0 .or. pdbref) then
749 read(inp,'(a)') pdbfile
750 if(me.eq.king.or..not.out1file)
751 & write (iout,'(2a)') 'PDB data will be read from file ',
752 & pdbfile(:ilen(pdbfile))
753 open(ipdbin,file=pdbfile,status='old',err=33)
755 33 write (iout,'(a)') 'Error opening PDB file.'
758 c print *,'Begin reading pdb data'
760 c print *,'Finished reading pdb data'
761 if(me.eq.king.or..not.out1file)
762 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
763 & ' nstart_sup=',nstart_sup
765 itype_pdb(i)=itype(i)
769 nct=nstart_sup+nsup-1
770 call contact(.false.,ncont_ref,icont_ref,co)
773 if(me.eq.king.or..not.out1file)
774 & write(iout,*)'Adding sidechains'
778 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
781 do while (fail.and.nsi.le.maxsi)
782 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
785 if(fail) write(iout,*)'Adding sidechain failed for res ',
786 & i,' after ',nsi,' trials'
791 if (indpdb.eq.0) then
792 C Read sequence if not taken from the pdb file.
794 c print *,'nres=',nres
795 if (iscode.gt.0) then
796 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
798 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
800 C Convert sequence to numeric code
802 itype(i)=rescode(i,sequence(i),iscode)
804 C Assign initial virtual bond lengths
810 vbld(i+nres)=dsc(iabs(itype(i)))
811 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
812 c write (iout,*) "i",i," itype",itype(i),
813 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
817 c print '(20i4)',(itype(i),i=1,nres)
820 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
822 if (itype(i).eq.ntyp1) then
826 else if (iabs(itype(i+1)).ne.20) then
828 else if (iabs(itype(i)).ne.20) then
835 if(me.eq.king.or..not.out1file)then
836 write (iout,*) "ITEL"
838 write (iout,*) i,itype(i),itel(i)
840 print *,'Call Read_Bridge.'
843 C 8/13/98 Set limits to generating the dihedral angles
848 read (inp,*) ndih_constr
849 if (ndih_constr.gt.0) then
851 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
852 if(me.eq.king.or..not.out1file)then
854 & 'There are',ndih_constr,' constraints on phi angles.'
856 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
860 phi0(i)=deg2rad*phi0(i)
861 drange(i)=deg2rad*drange(i)
863 if(me.eq.king.or..not.out1file)
864 & write (iout,*) 'FTORS',ftors
867 phibound(1,ii) = phi0(i)-drange(i)
868 phibound(2,ii) = phi0(i)+drange(i)
875 write (iout,'(a)') 'Boundaries in phi angle sampling:'
877 write (iout,'(a3,i5,2f10.1)')
878 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
884 cd print *,'NNT=',NNT,' NCT=',NCT
885 if (itype(1).eq.ntyp1) nnt=2
886 if (itype(nres).eq.ntyp1) nct=nct-1
888 if(me.eq.king.or..not.out1file)
889 & write (iout,'(a,i3)') 'nsup=',nsup
891 if (nsup.le.(nct-nnt+1)) then
892 do i=0,nct-nnt+1-nsup
893 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
899 & 'Error - sequences to be superposed do not match.'
902 do i=0,nsup-(nct-nnt+1)
903 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
905 nstart_sup=nstart_sup+i
911 & 'Error - sequences to be superposed do not match.'
914 if (nsup.eq.0) nsup=nct-nnt
915 if (nstart_sup.eq.0) nstart_sup=nnt
916 if (nstart_seq.eq.0) nstart_seq=nnt
917 if(me.eq.king.or..not.out1file)
918 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
919 & ' nstart_seq=',nstart_seq
921 c--- Zscore rms -------
922 if (nz_start.eq.0) nz_start=nnt
923 if (nz_end.eq.0 .and. nsup.gt.0) then
925 else if (nz_end.eq.0) then
928 if(me.eq.king.or..not.out1file)then
929 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
930 write (iout,*) 'IZ_SC=',iz_sc
932 c----------------------
935 if (.not.pdbref) then
936 call read_angles(inp,*38)
938 38 write (iout,'(a)') 'Error reading reference structure.'
940 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
941 stop 'Error reading reference structure'
945 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
955 call contact(.true.,ncont_ref,icont_ref,co)
957 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
959 if (constr_dist.gt.0) call read_dist_constr
960 write (iout,*) "After read_dist_constr nhpb",nhpb
961 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
963 if(me.eq.king.or..not.out1file)
964 & write (iout,*) 'Contact order:',co
966 if(me.eq.king.or..not.out1file)
967 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
970 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
972 if(me.eq.king.or..not.out1file)
973 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
974 & icont_ref(1,i),' ',
975 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
979 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
980 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
981 & modecalc.ne.10) then
982 C If input structure hasn't been supplied from the PDB file read or generate
984 if (iranconf.eq.0 .and. .not. extconf) then
985 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
986 & write (iout,'(a)') 'Initial geometry will be read in.'
988 read(inp,'(8f10.5)',end=36,err=36)
989 & ((c(l,k),l=1,3),k=1,nres),
990 & ((c(l,k+nres),l=1,3),k=nnt,nct)
991 write (iout,*) "Exit READ_CART"
992 write (iout,'(8f10.5)')
993 & ((c(l,k),l=1,3),k=1,nres),
994 & ((c(l,k+nres),l=1,3),k=nnt,nct)
995 call int_from_cart1(.true.)
996 write (iout,*) "Finish INT_TO_CART"
999 dc(j,i)=c(j,i+1)-c(j,i)
1000 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1004 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1006 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1007 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1013 call read_angles(inp,*36)
1016 36 write (iout,'(a)') 'Error reading angle file.'
1018 call mpi_finalize( MPI_COMM_WORLD,IERR )
1020 stop 'Error reading angle file.'
1022 else if (extconf) then
1023 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1024 & write (iout,'(a)') 'Extended chain initial geometry.'
1026 theta(i)=90d0*deg2rad
1029 phi(i)=180d0*deg2rad
1032 alph(i)=110d0*deg2rad
1035 omeg(i)=-120d0*deg2rad
1036 if (itype(i).le.0) omeg(i)=-omeg(i)
1039 if(me.eq.king.or..not.out1file)
1040 & write (iout,'(a)') 'Random-generated initial geometry.'
1044 if (me.eq.king .or. fg_rank.eq.0 .and. (
1045 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1049 call gen_rand_conf(itmp,*30)
1051 30 write (iout,*) 'Failed to generate random conformation',
1052 & ', itrial=',itrial
1053 write (*,*) 'Processor:',me,
1054 & ' Failed to generate random conformation',
1064 write (iout,'(a,i3,a)') 'Processor:',me,
1065 & ' error in generating random conformation.'
1066 write (*,'(a,i3,a)') 'Processor:',me,
1067 & ' error in generating random conformation.'
1070 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1075 & ' error in generating random conformation.'
1080 elseif (modecalc.eq.4) then
1081 read (inp,'(a)') intinname
1082 open (intin,file=intinname,status='old',err=333)
1083 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1084 & write (iout,'(a)') 'intinname',intinname
1085 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1087 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1089 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1091 stop 'Error opening angle file.'
1095 C Generate distance constraints, if the PDB structure is to be regularized.
1096 if (nthread.gt.0) then
1097 call read_threadbase
1100 if (me.eq.king .or. .not. out1file)
1102 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1103 write (iout,'(/a,i3,a)')
1104 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1105 write (iout,'(20i4)') (iss(i),i=1,ns)
1107 write(iout,*)"Running with dynamic disulfide-bond formation"
1109 write (iout,'(/a/)') 'Pre-formed links are:'
1115 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1116 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1122 if (ns.gt.0.and.dyn_ss) then
1126 forcon(i-nss)=forcon(i)
1133 dyn_ss_mask(iss(i))=.true.
1136 if (i2ndstr.gt.0) call secstrp2dihc
1137 c call geom_to_var(nvar,x)
1138 c call etotal(energia(0))
1139 c call enerprint(energia(0))
1140 c call briefout(0,etot)
1142 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1143 cd write (iout,'(a)') 'Variable list:'
1144 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1146 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1147 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1148 & 'Processor',myrank,': end reading molecular data.'
1152 c--------------------------------------------------------------------------
1153 logical function seq_comp(itypea,itypeb,length)
1155 integer length,itypea(length),itypeb(length)
1158 if (itypea(i).ne.itypeb(i)) then
1166 c-----------------------------------------------------------------------------
1167 subroutine read_bridge
1168 C Read information about disulfide bridges.
1169 implicit real*8 (a-h,o-z)
1170 include 'DIMENSIONS'
1174 include 'COMMON.IOUNITS'
1175 include 'COMMON.GEO'
1176 include 'COMMON.VAR'
1177 include 'COMMON.INTERACT'
1178 include 'COMMON.LOCAL'
1179 include 'COMMON.NAMES'
1180 include 'COMMON.CHAIN'
1181 include 'COMMON.FFIELD'
1182 include 'COMMON.SBRIDGE'
1183 include 'COMMON.HEADER'
1184 include 'COMMON.CONTROL'
1185 include 'COMMON.DBASE'
1186 include 'COMMON.THREAD'
1187 include 'COMMON.TIME1'
1188 include 'COMMON.SETUP'
1189 C Read bridging residues.
1190 read (inp,*) ns,(iss(i),i=1,ns)
1192 if(me.eq.king.or..not.out1file)
1193 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1194 C Check whether the specified bridging residues are cystines.
1196 if (itype(iss(i)).ne.1) then
1197 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1198 & 'Do you REALLY think that the residue ',
1199 & restyp(itype(iss(i))),i,
1200 & ' can form a disulfide bridge?!!!'
1201 write (*,'(2a,i3,a)')
1202 & 'Do you REALLY think that the residue ',
1203 & restyp(itype(iss(i))),i,
1204 & ' can form a disulfide bridge?!!!'
1206 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1211 C Read preformed bridges.
1213 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1215 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1218 C Check if the residues involved in bridges are in the specified list of
1219 C bridging residues.
1222 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1223 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1224 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1225 & ' contains residues present in other pairs.'
1226 write (*,'(a,i3,a)') 'Disulfide pair',i,
1227 & ' contains residues present in other pairs.'
1229 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1235 if (ihpb(i).eq.iss(j)) goto 10
1237 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1240 if (jhpb(i).eq.iss(j)) goto 20
1242 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1248 ihpb(i)=ihpb(i)+nres
1249 jhpb(i)=jhpb(i)+nres
1255 c----------------------------------------------------------------------------
1256 subroutine read_x(kanal,*)
1257 implicit real*8 (a-h,o-z)
1258 include 'DIMENSIONS'
1259 include 'COMMON.GEO'
1260 include 'COMMON.VAR'
1261 include 'COMMON.CHAIN'
1262 include 'COMMON.IOUNITS'
1263 include 'COMMON.CONTROL'
1264 include 'COMMON.LOCAL'
1265 include 'COMMON.INTERACT'
1266 c Read coordinates from input
1268 read(kanal,'(8f10.5)',end=10,err=10)
1269 & ((c(l,k),l=1,3),k=1,nres),
1270 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1273 c(j,2*nres)=c(j,nres)
1275 call int_from_cart1(.false.)
1278 dc(j,i)=c(j,i+1)-c(j,i)
1279 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1283 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1285 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1286 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1294 c----------------------------------------------------------------------------
1295 subroutine read_threadbase
1296 implicit real*8 (a-h,o-z)
1297 include 'DIMENSIONS'
1298 include 'COMMON.IOUNITS'
1299 include 'COMMON.GEO'
1300 include 'COMMON.VAR'
1301 include 'COMMON.INTERACT'
1302 include 'COMMON.LOCAL'
1303 include 'COMMON.NAMES'
1304 include 'COMMON.CHAIN'
1305 include 'COMMON.FFIELD'
1306 include 'COMMON.SBRIDGE'
1307 include 'COMMON.HEADER'
1308 include 'COMMON.CONTROL'
1309 include 'COMMON.DBASE'
1310 include 'COMMON.THREAD'
1311 include 'COMMON.TIME1'
1312 C Read pattern database for threading.
1313 read (icbase,*) nseq
1315 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1316 & nres_base(2,i),nres_base(3,i)
1317 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1319 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1320 c & nres_base(2,i),nres_base(3,i)
1321 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1325 if (weidis.eq.0.0D0) weidis=0.1D0
1334 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1335 write (iout,'(a,i5)') 'nexcl: ',nexcl
1336 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1339 c------------------------------------------------------------------------------
1340 subroutine setup_var
1341 implicit real*8 (a-h,o-z)
1342 include 'DIMENSIONS'
1343 include 'COMMON.IOUNITS'
1344 include 'COMMON.GEO'
1345 include 'COMMON.VAR'
1346 include 'COMMON.INTERACT'
1347 include 'COMMON.LOCAL'
1348 include 'COMMON.NAMES'
1349 include 'COMMON.CHAIN'
1350 include 'COMMON.FFIELD'
1351 include 'COMMON.SBRIDGE'
1352 include 'COMMON.HEADER'
1353 include 'COMMON.CONTROL'
1354 include 'COMMON.DBASE'
1355 include 'COMMON.THREAD'
1356 include 'COMMON.TIME1'
1357 C Set up variable list.
1363 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1365 ialph(i,1)=nvar+nside
1369 if (indphi.gt.0) then
1371 else if (indback.gt.0) then
1376 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1379 c----------------------------------------------------------------------------
1380 subroutine gen_dist_constr
1381 C Generate CA distance constraints.
1382 implicit real*8 (a-h,o-z)
1383 include 'DIMENSIONS'
1384 include 'COMMON.IOUNITS'
1385 include 'COMMON.GEO'
1386 include 'COMMON.VAR'
1387 include 'COMMON.INTERACT'
1388 include 'COMMON.LOCAL'
1389 include 'COMMON.NAMES'
1390 include 'COMMON.CHAIN'
1391 include 'COMMON.FFIELD'
1392 include 'COMMON.SBRIDGE'
1393 include 'COMMON.HEADER'
1394 include 'COMMON.CONTROL'
1395 include 'COMMON.DBASE'
1396 include 'COMMON.THREAD'
1397 include 'COMMON.TIME1'
1398 dimension itype_pdb(maxres)
1399 common /pizda/ itype_pdb
1401 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1402 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1403 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1405 do i=nstart_sup,nstart_sup+nsup-1
1406 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1407 cd & ' seq_pdb', restyp(itype_pdb(i))
1408 do j=i+2,nstart_sup+nsup-1
1410 ihpb(nhpb)=i+nstart_seq-nstart_sup
1411 jhpb(nhpb)=j+nstart_seq-nstart_sup
1413 dhpb(nhpb)=dist(i,j)
1416 cd write (iout,'(a)') 'Distance constraints:'
1421 cd if (ii.gt.nres) then
1426 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1427 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1428 cd & dhpb(i),forcon(i)
1432 c----------------------------------------------------------------------------
1434 implicit real*8 (a-h,o-z)
1435 include 'DIMENSIONS'
1436 include 'COMMON.MAP'
1437 include 'COMMON.IOUNITS'
1438 character*3 angid(4) /'THE','PHI','ALP','OME'/
1439 character*80 mapcard,ucase
1441 read (inp,'(a)') mapcard
1442 mapcard=ucase(mapcard)
1443 if (index(mapcard,'PHI').gt.0) then
1445 else if (index(mapcard,'THE').gt.0) then
1447 else if (index(mapcard,'ALP').gt.0) then
1449 else if (index(mapcard,'OME').gt.0) then
1452 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1453 stop 'Error - illegal variable spec in MAP card.'
1455 call readi (mapcard,'RES1',res1(imap),0)
1456 call readi (mapcard,'RES2',res2(imap),0)
1457 if (res1(imap).eq.0) then
1458 res1(imap)=res2(imap)
1459 else if (res2(imap).eq.0) then
1460 res2(imap)=res1(imap)
1462 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1464 & 'Error - illegal definition of variable group in MAP.'
1465 stop 'Error - illegal definition of variable group in MAP.'
1467 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1468 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1469 call readi(mapcard,'NSTEP',nstep(imap),0)
1470 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1472 & 'Illegal boundary and/or step size specification in MAP.'
1473 stop 'Illegal boundary and/or step size specification in MAP.'
1478 c----------------------------------------------------------------------------
1480 implicit real*8 (a-h,o-z)
1481 include 'DIMENSIONS'
1482 include 'COMMON.IOUNITS'
1483 include 'COMMON.GEO'
1484 include 'COMMON.CSA'
1485 include 'COMMON.BANK'
1486 include 'COMMON.CONTROL'
1488 character*620 mcmcard
1489 call card_concat(mcmcard)
1491 call readi(mcmcard,'NCONF',nconf,50)
1492 call readi(mcmcard,'NADD',nadd,0)
1493 call readi(mcmcard,'JSTART',jstart,1)
1494 call readi(mcmcard,'JEND',jend,1)
1495 call readi(mcmcard,'NSTMAX',nstmax,500000)
1496 call readi(mcmcard,'N0',n0,1)
1497 call readi(mcmcard,'N1',n1,6)
1498 call readi(mcmcard,'N2',n2,4)
1499 call readi(mcmcard,'N3',n3,0)
1500 call readi(mcmcard,'N4',n4,0)
1501 call readi(mcmcard,'N5',n5,0)
1502 call readi(mcmcard,'N6',n6,10)
1503 call readi(mcmcard,'N7',n7,0)
1504 call readi(mcmcard,'N8',n8,0)
1505 call readi(mcmcard,'N9',n9,0)
1506 call readi(mcmcard,'N14',n14,0)
1507 call readi(mcmcard,'N15',n15,0)
1508 call readi(mcmcard,'N16',n16,0)
1509 call readi(mcmcard,'N17',n17,0)
1510 call readi(mcmcard,'N18',n18,0)
1512 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1514 call readi(mcmcard,'NDIFF',ndiff,2)
1515 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1516 call readi(mcmcard,'IS1',is1,1)
1517 call readi(mcmcard,'IS2',is2,8)
1518 call readi(mcmcard,'NRAN0',nran0,4)
1519 call readi(mcmcard,'NRAN1',nran1,2)
1520 call readi(mcmcard,'IRR',irr,1)
1521 call readi(mcmcard,'NSEED',nseed,20)
1522 call readi(mcmcard,'NTOTAL',ntotal,10000)
1523 call reada(mcmcard,'CUT1',cut1,2.0d0)
1524 call reada(mcmcard,'CUT2',cut2,5.0d0)
1525 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1526 call readi(mcmcard,'ICMAX',icmax,3)
1527 call readi(mcmcard,'IRESTART',irestart,0)
1528 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1531 call reada(mcmcard,'DELE',dele,20.0d0)
1532 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1533 call readi(mcmcard,'IREF',iref,0)
1534 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1535 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1536 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1537 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1538 write (iout,*) "NCONF_IN",nconf_in
1541 c----------------------------------------------------------------------------
1542 cfmc subroutine mcmfread
1543 cfmc implicit real*8 (a-h,o-z)
1544 cfmc include 'DIMENSIONS'
1545 cfmc include 'COMMON.MCMF'
1546 cfmc include 'COMMON.IOUNITS'
1547 cfmc include 'COMMON.GEO'
1548 cfmc character*80 ucase
1549 cfmc character*620 mcmcard
1550 cfmc call card_concat(mcmcard)
1552 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1553 cfmc write(iout,*)'MAXRANT=',maxrant
1554 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1555 cfmc write(iout,*)'MAXFAM=',maxfam
1556 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1557 cfmc write(iout,*)'NNET1=',nnet1
1558 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1559 cfmc write(iout,*)'NNET2=',nnet2
1560 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1561 cfmc write(iout,*)'NNET3=',nnet3
1562 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1563 cfmc write(iout,*)'ILASTT=',ilastt
1564 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1565 cfmc write(iout,*)'MAXSTR=',maxstr
1566 cfmc maxstr_f=maxstr/maxfam
1567 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1568 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1569 cfmc write(iout,*)'NMCMF=',nmcmf
1570 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1571 cfmc write(iout,*)'IFOCUS=',ifocus
1572 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1573 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1574 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1575 cfmc write(iout,*)'INTPRT=',intprt
1576 cfmc call readi(mcmcard,'IPRT',iprt,100)
1577 cfmc write(iout,*)'IPRT=',iprt
1578 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1579 cfmc write(iout,*)'IMAXTR=',imaxtr
1580 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1581 cfmc write(iout,*)'MAXEVEN=',maxeven
1582 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1583 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1584 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1585 cfmc write(iout,*)'INIMIN=',inimin
1586 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1587 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1588 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1589 cfmc write(iout,*)'NTHREAD=',nthread
1590 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1591 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1592 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1593 cfmc write(iout,*)'MAXPERT=',maxpert
1594 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1595 cfmc write(iout,*)'IRMSD=',irmsd
1596 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1597 cfmc write(iout,*)'DENEMIN=',denemin
1598 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1599 cfmc write(iout,*)'RCUT1S=',rcut1s
1600 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1601 cfmc write(iout,*)'RCUT1E=',rcut1e
1602 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1603 cfmc write(iout,*)'RCUT2S=',rcut2s
1604 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1605 cfmc write(iout,*)'RCUT2E=',rcut2e
1606 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1607 cfmc write(iout,*)'DPERT1=',d_pert1
1608 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1609 cfmc write(iout,*)'DPERT1A=',d_pert1a
1610 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1611 cfmc write(iout,*)'DPERT2=',d_pert2
1612 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1613 cfmc write(iout,*)'DPERT2A=',d_pert2a
1614 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1615 cfmc write(iout,*)'DPERT2B=',d_pert2b
1616 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1617 cfmc write(iout,*)'DPERT2C=',d_pert2c
1618 cfmc d_pert1=deg2rad*d_pert1
1619 cfmc d_pert1a=deg2rad*d_pert1a
1620 cfmc d_pert2=deg2rad*d_pert2
1621 cfmc d_pert2a=deg2rad*d_pert2a
1622 cfmc d_pert2b=deg2rad*d_pert2b
1623 cfmc d_pert2c=deg2rad*d_pert2c
1624 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1625 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1626 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1627 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1628 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1629 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1630 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1631 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1632 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1633 cfmc write(iout,*)'RCUTINI=',rcutini
1634 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1635 cfmc write(iout,*)'GRAT=',grat
1636 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1637 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1641 c----------------------------------------------------------------------------
1643 implicit real*8 (a-h,o-z)
1644 include 'DIMENSIONS'
1645 include 'COMMON.MCM'
1646 include 'COMMON.MCE'
1647 include 'COMMON.IOUNITS'
1649 character*320 mcmcard
1650 call card_concat(mcmcard)
1651 call readi(mcmcard,'MAXACC',maxacc,100)
1652 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1653 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1654 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1655 call readi(mcmcard,'MAXREPM',maxrepm,200)
1656 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1657 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1658 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1659 call reada(mcmcard,'E_UP',e_up,5.0D0)
1660 call reada(mcmcard,'DELTE',delte,0.1D0)
1661 call readi(mcmcard,'NSWEEP',nsweep,5)
1662 call readi(mcmcard,'NSTEPH',nsteph,0)
1663 call readi(mcmcard,'NSTEPC',nstepc,0)
1664 call reada(mcmcard,'TMIN',tmin,298.0D0)
1665 call reada(mcmcard,'TMAX',tmax,298.0D0)
1666 call readi(mcmcard,'NWINDOW',nwindow,0)
1667 call readi(mcmcard,'PRINT_MC',print_mc,0)
1668 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1669 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1670 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1671 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1672 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1673 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1674 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1675 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1676 if (nwindow.gt.0) then
1677 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1679 winlen(i)=winend(i)-winstart(i)+1
1682 if (tmax.lt.tmin) tmax=tmin
1683 if (tmax.eq.tmin) then
1687 if (nstepc.gt.0 .and. nsteph.gt.0) then
1688 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1689 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1691 C Probabilities of different move types
1692 sumpro_type(0)=0.0D0
1693 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1694 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1695 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1696 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1697 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1698 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1699 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1701 print *,'i',i,' sumprotype',sumpro_type(i)
1702 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1703 print *,'i',i,' sumprotype',sumpro_type(i)
1707 c----------------------------------------------------------------------------
1708 subroutine read_minim
1709 implicit real*8 (a-h,o-z)
1710 include 'DIMENSIONS'
1711 include 'COMMON.MINIM'
1712 include 'COMMON.IOUNITS'
1714 character*320 minimcard
1715 call card_concat(minimcard)
1716 call readi(minimcard,'MAXMIN',maxmin,2000)
1717 call readi(minimcard,'MAXFUN',maxfun,5000)
1718 call readi(minimcard,'MINMIN',minmin,maxmin)
1719 call readi(minimcard,'MINFUN',minfun,maxmin)
1720 call reada(minimcard,'TOLF',tolf,1.0D-2)
1721 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1722 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1723 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1724 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1725 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1726 & 'Options in energy minimization:'
1727 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1728 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1729 & 'MinMin:',MinMin,' MinFun:',MinFun,
1730 & ' TolF:',TolF,' RTolF:',RTolF
1733 c----------------------------------------------------------------------------
1734 subroutine read_angles(kanal,*)
1735 implicit real*8 (a-h,o-z)
1736 include 'DIMENSIONS'
1737 include 'COMMON.GEO'
1738 include 'COMMON.VAR'
1739 include 'COMMON.CHAIN'
1740 include 'COMMON.IOUNITS'
1741 include 'COMMON.CONTROL'
1742 c Read angles from input
1744 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1745 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1746 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1747 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1750 c 9/7/01 avoid 180 deg valence angle
1751 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1753 theta(i)=deg2rad*theta(i)
1754 phi(i)=deg2rad*phi(i)
1755 alph(i)=deg2rad*alph(i)
1756 omeg(i)=deg2rad*omeg(i)
1761 c----------------------------------------------------------------------------
1762 subroutine reada(rekord,lancuch,wartosc,default)
1764 character*(*) rekord,lancuch
1765 double precision wartosc,default
1768 iread=index(rekord,lancuch)
1769 if (iread.eq.0) then
1773 iread=iread+ilen(lancuch)+1
1774 read (rekord(iread:),*,err=10,end=10) wartosc
1779 c----------------------------------------------------------------------------
1780 subroutine readi(rekord,lancuch,wartosc,default)
1782 character*(*) rekord,lancuch
1783 integer wartosc,default
1786 iread=index(rekord,lancuch)
1787 if (iread.eq.0) then
1791 iread=iread+ilen(lancuch)+1
1792 read (rekord(iread:),*,err=10,end=10) wartosc
1797 c----------------------------------------------------------------------------
1798 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1801 integer tablica(dim),default
1802 character*(*) rekord,lancuch
1809 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1810 if (iread.eq.0) return
1811 iread=iread+ilen(lancuch)+1
1812 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1815 c----------------------------------------------------------------------------
1816 subroutine multreada(rekord,lancuch,tablica,dim,default)
1819 double precision tablica(dim),default
1820 character*(*) rekord,lancuch
1827 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1828 if (iread.eq.0) return
1829 iread=iread+ilen(lancuch)+1
1830 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1833 c----------------------------------------------------------------------------
1834 subroutine openunits
1835 implicit real*8 (a-h,o-z)
1836 include 'DIMENSIONS'
1839 character*16 form,nodename
1842 include 'COMMON.SETUP'
1843 include 'COMMON.IOUNITS'
1845 include 'COMMON.CONTROL'
1846 integer lenpre,lenpot,ilen,lentmp
1848 character*3 out1file_text,ucase
1851 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1852 call getenv_loc("PREFIX",prefix)
1854 call getenv_loc("POT",pot)
1855 call getenv_loc("DIRTMP",tmpdir)
1856 call getenv_loc("CURDIR",curdir)
1857 call getenv_loc("OUT1FILE",out1file_text)
1858 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1859 out1file_text=ucase(out1file_text)
1860 if (out1file_text(1:1).eq."Y") then
1863 out1file=fg_rank.gt.0
1868 if (lentmp.gt.0) then
1869 write (*,'(80(1h!))')
1870 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1871 write (*,'(80(1h!))')
1872 write (*,*)"All output files will be on node /tmp directory."
1874 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1875 if (me.eq.king) then
1876 write (*,*) "The master node is ",nodename
1877 else if (fg_rank.eq.0) then
1878 write (*,*) "I am the CG slave node ",nodename
1880 write (*,*) "I am the FG slave node ",nodename
1883 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1884 lenpre = lentmp+lenpre+1
1886 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1887 C Get the names and open the input files
1888 #if defined(WINIFL) || defined(WINPGI)
1889 open(1,file=pref_orig(:ilen(pref_orig))//
1890 & '.inp',status='old',readonly,shared)
1891 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1892 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1893 C Get parameter filenames and open the parameter files.
1894 call getenv_loc('BONDPAR',bondname)
1895 open (ibond,file=bondname,status='old',readonly,shared)
1896 call getenv_loc('THETPAR',thetname)
1897 open (ithep,file=thetname,status='old',readonly,shared)
1898 call getenv_loc('ROTPAR',rotname)
1899 open (irotam,file=rotname,status='old',readonly,shared)
1900 call getenv_loc('TORPAR',torname)
1901 open (itorp,file=torname,status='old',readonly,shared)
1902 call getenv_loc('TORDPAR',tordname)
1903 open (itordp,file=tordname,status='old',readonly,shared)
1904 call getenv_loc('FOURIER',fouriername)
1905 open (ifourier,file=fouriername,status='old',readonly,shared)
1906 call getenv_loc('ELEPAR',elename)
1907 open (ielep,file=elename,status='old',readonly,shared)
1908 call getenv_loc('SIDEPAR',sidename)
1909 open (isidep,file=sidename,status='old',readonly,shared)
1910 #elif (defined CRAY) || (defined AIX)
1911 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1913 c print *,"Processor",myrank," opened file 1"
1914 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1915 c print *,"Processor",myrank," opened file 9"
1916 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1917 C Get parameter filenames and open the parameter files.
1918 call getenv_loc('BONDPAR',bondname)
1919 open (ibond,file=bondname,status='old',action='read')
1920 c print *,"Processor",myrank," opened file IBOND"
1921 call getenv_loc('THETPAR',thetname)
1922 open (ithep,file=thetname,status='old',action='read')
1923 c print *,"Processor",myrank," opened file ITHEP"
1924 call getenv_loc('ROTPAR',rotname)
1925 open (irotam,file=rotname,status='old',action='read')
1926 c print *,"Processor",myrank," opened file IROTAM"
1927 call getenv_loc('TORPAR',torname)
1928 open (itorp,file=torname,status='old',action='read')
1929 c print *,"Processor",myrank," opened file ITORP"
1930 call getenv_loc('TORDPAR',tordname)
1931 open (itordp,file=tordname,status='old',action='read')
1932 c print *,"Processor",myrank," opened file ITORDP"
1933 call getenv_loc('SCCORPAR',sccorname)
1934 open (isccor,file=sccorname,status='old',action='read')
1935 c print *,"Processor",myrank," opened file ISCCOR"
1936 call getenv_loc('FOURIER',fouriername)
1937 open (ifourier,file=fouriername,status='old',action='read')
1938 c print *,"Processor",myrank," opened file IFOURIER"
1939 call getenv_loc('ELEPAR',elename)
1940 open (ielep,file=elename,status='old',action='read')
1941 c print *,"Processor",myrank," opened file IELEP"
1942 call getenv_loc('SIDEPAR',sidename)
1943 open (isidep,file=sidename,status='old',action='read')
1944 c print *,"Processor",myrank," opened file ISIDEP"
1945 c print *,"Processor",myrank," opened parameter files"
1947 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1948 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1949 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1950 C Get parameter filenames and open the parameter files.
1951 call getenv_loc('BONDPAR',bondname)
1952 open (ibond,file=bondname,status='old')
1953 call getenv_loc('THETPAR',thetname)
1954 open (ithep,file=thetname,status='old')
1955 call getenv_loc('ROTPAR',rotname)
1956 open (irotam,file=rotname,status='old')
1957 call getenv_loc('TORPAR',torname)
1958 open (itorp,file=torname,status='old')
1959 call getenv_loc('TORDPAR',tordname)
1960 open (itordp,file=tordname,status='old')
1961 call getenv_loc('SCCORPAR',sccorname)
1962 open (isccor,file=sccorname,status='old')
1963 call getenv_loc('FOURIER',fouriername)
1964 open (ifourier,file=fouriername,status='old')
1965 call getenv_loc('ELEPAR',elename)
1966 open (ielep,file=elename,status='old')
1967 call getenv_loc('SIDEPAR',sidename)
1968 open (isidep,file=sidename,status='old')
1970 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1972 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1973 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1974 C Get parameter filenames and open the parameter files.
1975 call getenv_loc('BONDPAR',bondname)
1976 open (ibond,file=bondname,status='old',readonly)
1977 call getenv_loc('THETPAR',thetname)
1978 open (ithep,file=thetname,status='old',readonly)
1979 call getenv_loc('ROTPAR',rotname)
1980 open (irotam,file=rotname,status='old',readonly)
1981 call getenv_loc('TORPAR',torname)
1982 open (itorp,file=torname,status='old',readonly)
1983 call getenv_loc('TORDPAR',tordname)
1984 open (itordp,file=tordname,status='old',readonly)
1985 call getenv_loc('SCCORPAR',sccorname)
1986 open (isccor,file=sccorname,status='old',readonly)
1988 call getenv_loc('THETPARPDB',thetname_pdb)
1989 print *,"thetname_pdb ",thetname_pdb
1990 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1991 print *,ithep_pdb," opened"
1993 call getenv_loc('FOURIER',fouriername)
1994 open (ifourier,file=fouriername,status='old',readonly)
1995 call getenv_loc('ELEPAR',elename)
1996 open (ielep,file=elename,status='old',readonly)
1997 call getenv_loc('SIDEPAR',sidename)
1998 open (isidep,file=sidename,status='old',readonly)
1999 call getenv_loc('LIPTRANPAR',liptranname)
2000 open (iliptranpar,file=liptranname,status='old',action='read')
2002 call getenv_loc('ROTPARPDB',rotname_pdb)
2003 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2008 C 8/9/01 In the newest version SCp interaction constants are read from a file
2009 C Use -DOLDSCP to use hard-coded constants instead.
2011 call getenv_loc('SCPPAR',scpname)
2012 #if defined(WINIFL) || defined(WINPGI)
2013 open (iscpp,file=scpname,status='old',readonly,shared)
2014 #elif (defined CRAY) || (defined AIX)
2015 open (iscpp,file=scpname,status='old',action='read')
2017 open (iscpp,file=scpname,status='old')
2019 open (iscpp,file=scpname,status='old',readonly)
2022 call getenv_loc('PATTERN',patname)
2023 #if defined(WINIFL) || defined(WINPGI)
2024 open (icbase,file=patname,status='old',readonly,shared)
2025 #elif (defined CRAY) || (defined AIX)
2026 open (icbase,file=patname,status='old',action='read')
2028 open (icbase,file=patname,status='old')
2030 open (icbase,file=patname,status='old',readonly)
2033 C Open output file only for CG processes
2034 c print *,"Processor",myrank," fg_rank",fg_rank
2035 if (fg_rank.eq.0) then
2037 if (nodes.eq.1) then
2040 npos = dlog10(dfloat(nodes-1))+1
2042 if (npos.lt.3) npos=3
2043 write (liczba,'(i1)') npos
2044 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2046 write (liczba,form) me
2047 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2048 & liczba(:ilen(liczba))
2049 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2051 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2053 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2054 & liczba(:ilen(liczba))//'.mol2'
2055 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2056 & liczba(:ilen(liczba))//'.stat'
2058 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2059 & //liczba(:ilen(liczba))//'.stat')
2060 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2063 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2064 & liczba(:ilen(liczba))//'.const'
2069 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2070 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2071 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2072 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2073 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2075 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2077 rest2name=prefix(:ilen(prefix))//'.rst'
2079 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2082 #if defined(AIX) || defined(PGI)
2083 if (me.eq.king .or. .not. out1file)
2084 & open(iout,file=outname,status='unknown')
2086 if (fg_rank.gt.0) then
2087 write (liczba,'(i3.3)') myrank/nfgtasks
2088 write (ll,'(bz,i3.3)') fg_rank
2089 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2094 open(igeom,file=intname,status='unknown',position='append')
2095 open(ipdb,file=pdbname,status='unknown')
2096 open(imol2,file=mol2name,status='unknown')
2097 open(istat,file=statname,status='unknown',position='append')
2099 c1out open(iout,file=outname,status='unknown')
2102 if (me.eq.king .or. .not.out1file)
2103 & open(iout,file=outname,status='unknown')
2105 if (fg_rank.gt.0) then
2106 write (liczba,'(i3.3)') myrank/nfgtasks
2107 write (ll,'(bz,i3.3)') fg_rank
2108 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2113 open(igeom,file=intname,status='unknown',access='append')
2114 open(ipdb,file=pdbname,status='unknown')
2115 open(imol2,file=mol2name,status='unknown')
2116 open(istat,file=statname,status='unknown',access='append')
2118 c1out open(iout,file=outname,status='unknown')
2121 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2122 csa_seed=prefix(:lenpre)//'.CSA.seed'
2123 csa_history=prefix(:lenpre)//'.CSA.history'
2124 csa_bank=prefix(:lenpre)//'.CSA.bank'
2125 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2126 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2127 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2128 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2129 csa_int=prefix(:lenpre)//'.int'
2130 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2131 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2132 csa_in=prefix(:lenpre)//'.CSA.in'
2133 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2136 write (iout,'(80(1h-))')
2137 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2138 write (iout,'(80(1h-))')
2139 write (iout,*) "Input file : ",
2140 & pref_orig(:ilen(pref_orig))//'.inp'
2141 write (iout,*) "Output file : ",
2142 & outname(:ilen(outname))
2144 write (iout,*) "Sidechain potential file : ",
2145 & sidename(:ilen(sidename))
2147 write (iout,*) "SCp potential file : ",
2148 & scpname(:ilen(scpname))
2150 write (iout,*) "Electrostatic potential file : ",
2151 & elename(:ilen(elename))
2152 write (iout,*) "Cumulant coefficient file : ",
2153 & fouriername(:ilen(fouriername))
2154 write (iout,*) "Torsional parameter file : ",
2155 & torname(:ilen(torname))
2156 write (iout,*) "Double torsional parameter file : ",
2157 & tordname(:ilen(tordname))
2158 write (iout,*) "SCCOR parameter file : ",
2159 & sccorname(:ilen(sccorname))
2160 write (iout,*) "Bond & inertia constant file : ",
2161 & bondname(:ilen(bondname))
2162 write (iout,*) "Bending parameter file : ",
2163 & thetname(:ilen(thetname))
2164 write (iout,*) "Rotamer parameter file : ",
2165 & rotname(:ilen(rotname))
2166 write (iout,*) "Thetpdb parameter file : ",
2167 & thetname_pdb(:ilen(thetname_pdb))
2168 write (iout,*) "Threading database : ",
2169 & patname(:ilen(patname))
2171 &write (iout,*)" DIRTMP : ",
2173 write (iout,'(80(1h-))')
2177 c----------------------------------------------------------------------------
2178 subroutine card_concat(card)
2179 implicit real*8 (a-h,o-z)
2180 include 'DIMENSIONS'
2181 include 'COMMON.IOUNITS'
2183 character*80 karta,ucase
2185 read (inp,'(a)') karta
2188 do while (karta(80:80).eq.'&')
2189 card=card(:ilen(card)+1)//karta(:79)
2190 read (inp,'(a)') karta
2193 card=card(:ilen(card)+1)//karta
2196 c----------------------------------------------------------------------------------
2198 implicit real*8 (a-h,o-z)
2199 include 'DIMENSIONS'
2200 include 'COMMON.CHAIN'
2201 include 'COMMON.IOUNITS'
2203 open(irest2,file=rest2name,status='unknown')
2204 read(irest2,*) totT,EK,potE,totE,t_bath
2207 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2210 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2213 read (irest2,*) iset
2218 c---------------------------------------------------------------------------------
2219 subroutine read_fragments
2220 implicit real*8 (a-h,o-z)
2221 include 'DIMENSIONS'
2225 include 'COMMON.SETUP'
2226 include 'COMMON.CHAIN'
2227 include 'COMMON.IOUNITS'
2229 include 'COMMON.CONTROL'
2230 read(inp,*) nset,nfrag,npair,nfrag_back
2231 if(me.eq.king.or..not.out1file)
2232 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2233 & " nfrag_back",nfrag_back
2235 read(inp,*) mset(iset)
2237 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2239 if(me.eq.king.or..not.out1file)
2240 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2241 & ifrag(2,i,iset), qinfrag(i,iset)
2244 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2246 if(me.eq.king.or..not.out1file)
2247 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2248 & ipair(2,i,iset), qinpair(i,iset)
2251 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2252 & wfrag_back(3,i,iset),
2253 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2254 if(me.eq.king.or..not.out1file)
2255 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2256 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2261 C---------------------------------------------------------------------------
2262 subroutine read_afminp
2263 implicit real*8 (a-h,o-z)
2264 include 'DIMENSIONS'
2268 include 'COMMON.SETUP'
2269 include 'COMMON.CONTROL'
2270 include 'COMMON.CHAIN'
2271 include 'COMMON.IOUNITS'
2272 include 'COMMON.SBRIDGE'
2273 character*320 afmcard
2275 call card_concat(afmcard)
2276 call readi(afmcard,"BEG",afmbeg,0)
2277 call readi(afmcard,"END",afmend,0)
2278 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2279 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2280 print *,'FORCE=' ,forceAFMconst
2281 CCCC NOW PROPERTIES FOR AFM
2284 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2286 distafminit=dsqrt(distafminit)
2287 print *,'initdist',distafminit
2291 c-------------------------------------------------------------------------------
2292 subroutine read_dist_constr
2293 implicit real*8 (a-h,o-z)
2294 include 'DIMENSIONS'
2298 include 'COMMON.SETUP'
2299 include 'COMMON.CONTROL'
2300 include 'COMMON.CHAIN'
2301 include 'COMMON.IOUNITS'
2302 include 'COMMON.SBRIDGE'
2303 integer ifrag_(2,100),ipair_(2,100)
2304 double precision wfrag_(100),wpair_(100)
2305 character*500 controlcard
2306 c write (iout,*) "Calling read_dist_constr"
2307 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2309 call card_concat(controlcard)
2310 call readi(controlcard,"NFRAG",nfrag_,0)
2311 call readi(controlcard,"NPAIR",npair_,0)
2312 call readi(controlcard,"NDIST",ndist_,0)
2313 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2314 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2315 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2316 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2317 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2318 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2319 c write (iout,*) "IFRAG"
2321 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2323 c write (iout,*) "IPAIR"
2325 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2329 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2330 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2331 & ifrag_(2,i)=nstart_sup+nsup-1
2332 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2334 if (wfrag_(i).gt.0.0d0) then
2335 do j=ifrag_(1,i),ifrag_(2,i)-1
2336 do k=j+1,ifrag_(2,i)
2337 c write (iout,*) "j",j," k",k
2339 if (constr_dist.eq.1) then
2344 forcon(nhpb)=wfrag_(i)
2345 else if (constr_dist.eq.2) then
2346 if (ddjk.le.dist_cut) then
2351 forcon(nhpb)=wfrag_(i)
2358 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2361 if (.not.out1file .or. me.eq.king)
2362 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2363 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2365 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2366 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2373 if (wpair_(i).gt.0.0d0) then
2381 do j=ifrag_(1,ii),ifrag_(2,ii)
2382 do k=ifrag_(1,jj),ifrag_(2,jj)
2386 forcon(nhpb)=wpair_(i)
2387 dhpb(nhpb)=dist(j,k)
2389 if (.not.out1file .or. me.eq.king)
2390 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2391 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2393 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2394 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2401 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2402 if (forcon(nhpb+1).gt.0.0d0) then
2404 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2406 if (.not.out1file .or. me.eq.king)
2407 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2408 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2410 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2411 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2418 c-------------------------------------------------------------------------------
2420 subroutine flush(iu)
2425 subroutine flush(iu)
2430 c------------------------------------------------------------------------------
2431 subroutine copy_to_tmp(source)
2432 include "DIMENSIONS"
2433 include "COMMON.IOUNITS"
2434 character*(*) source
2435 character* 256 tmpfile
2439 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2440 inquire(file=tmpfile,exist=ex)
2442 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2443 & " to temporary directory..."
2444 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2445 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2449 c------------------------------------------------------------------------------
2450 subroutine move_from_tmp(source)
2451 include "DIMENSIONS"
2452 include "COMMON.IOUNITS"
2453 character*(*) source
2456 write (*,*) "Moving ",source(:ilen(source)),
2457 & " from temporary directory to working directory"
2458 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2459 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2462 c------------------------------------------------------------------------------
2463 subroutine random_init(seed)
2465 C Initialize random number generator
2467 implicit real*8 (a-h,o-z)
2468 include 'DIMENSIONS'
2471 logical OKRandom, prng_restart
2473 integer iseed_array(4)
2475 include 'COMMON.IOUNITS'
2476 include 'COMMON.TIME1'
2477 include 'COMMON.THREAD'
2478 include 'COMMON.SBRIDGE'
2479 include 'COMMON.CONTROL'
2480 include 'COMMON.MCM'
2481 include 'COMMON.MAP'
2482 include 'COMMON.HEADER'
2483 include 'COMMON.CSA'
2484 include 'COMMON.CHAIN'
2485 include 'COMMON.MUCA'
2487 include 'COMMON.FFIELD'
2488 include 'COMMON.SETUP'
2489 iseed=-dint(dabs(seed))
2490 if (iseed.eq.0) then
2491 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2492 & 'Random seed undefined. The program will stop.'
2493 write (*,'(/80(1h*)/20x,a/80(1h*))')
2494 & 'Random seed undefined. The program will stop.'
2496 call mpi_finalize(mpi_comm_world,ierr)
2498 stop 'Bad random seed.'
2501 if (fg_rank.eq.0) then
2505 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2506 OKRandom = prng_restart(me,iseed)
2509 tmp=65536.0d0**(4-i)
2510 iseed_array(i) = dint(seed/tmp)
2511 seed=seed-iseed_array(i)*tmp
2514 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2515 & (iseed_array(i),i=1,4)
2516 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2517 & (iseed_array(i),i=1,4)
2518 OKRandom = prng_restart(me,iseed_array)
2521 c r1 = prng_next(me)
2522 r1=ran_number(0.0D0,1.0D0)
2524 & write (iout,*) 'ran_num',r1
2525 if (r1.lt.0.0d0) OKRandom=.false.
2527 if (.not.OKRandom) then
2528 write (iout,*) 'PRNG IS NOT WORKING!!!'
2529 print *,'PRNG IS NOT WORKING!!!'
2532 call mpi_abort(mpi_comm_world,error_msg,ierr)
2535 write (iout,*) 'too many processors for parallel prng'
2536 write (*,*) 'too many processors for parallel prng'
2544 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)