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 include 'COMMON.SHIELD'
85 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
86 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
88 character*320 controlcard
93 read (INP,'(a)') titel
94 call card_concat(controlcard)
95 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
96 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
97 call reada(controlcard,'SEED',seed,0.0D0)
98 call random_init(seed)
99 C Set up the time limit (caution! The time must be input in minutes!)
100 read_cart=index(controlcard,'READ_CART').gt.0
101 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
102 C this variable with_theta_constr is the variable which allow to read and execute the
103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
104 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
105 write (iout,*) "with_theta_constr ",with_theta_constr
106 call readi(controlcard,'SYM',symetr,1)
107 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
108 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
109 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
110 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
111 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
112 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
113 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
114 call reada(controlcard,'DRMS',drms,0.1D0)
115 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
117 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
118 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
119 write (iout,'(a,f10.1)')'DRMS = ',drms
120 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
121 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
123 call readi(controlcard,'NZ_START',nz_start,0)
124 call readi(controlcard,'NZ_END',nz_end,0)
125 call readi(controlcard,'IZ_SC',iz_sc,0)
127 safety = 60.0d0*safety
130 call reada(controlcard,"T_BATH",t_bath,300.0d0)
131 minim=(index(controlcard,'MINIMIZE').gt.0)
132 dccart=(index(controlcard,'CART').gt.0)
133 overlapsc=(index(controlcard,'OVERLAP').gt.0)
134 overlapsc=.not.overlapsc
135 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
136 searchsc=.not.searchsc
137 sideadd=(index(controlcard,'SIDEADD').gt.0)
138 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
139 outpdb=(index(controlcard,'PDBOUT').gt.0)
140 outmol2=(index(controlcard,'MOL2OUT').gt.0)
141 pdbref=(index(controlcard,'PDBREF').gt.0)
142 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
143 indpdb=index(controlcard,'PDBSTART')
144 extconf=(index(controlcard,'EXTCONF').gt.0)
145 AFMlog=(index(controlcard,'AFM'))
146 selfguide=(index(controlcard,'SELFGUIDE'))
147 print *,'AFMlog',AFMlog,selfguide,"KUPA"
148 call readi(controlcard,'TUBEMOD',tubelog,0)
149 write (iout,*) TUBElog,"TUBEMODE"
150 call readi(controlcard,'GENCONSTR',genconstr,0)
151 C write (iout,*) TUBElog,"TUBEMODE"
152 call readi(controlcard,'IPRINT',iprint,0)
153 C SHIELD keyword sets if the shielding effect of side-chains is used
154 C 0 denots no shielding is used all peptide are equally despite the
155 C solvent accesible area
156 C 1 the newly introduced function
157 C 2 reseved for further possible developement
158 call readi(controlcard,'SHIELD',shield_mode,0)
159 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
160 write(iout,*) "shield_mode",shield_mode
162 call readi(controlcard,'TORMODE',tor_mode,0)
163 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
164 write(iout,*) "torsional and valence angle mode",tor_mode
165 call readi(controlcard,'MAXGEN',maxgen,10000)
166 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
167 call readi(controlcard,"KDIAG",kdiag,0)
168 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
169 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
170 & write (iout,*) "RESCALE_MODE",rescale_mode
171 split_ene=index(controlcard,'SPLIT_ENE').gt.0
172 if (index(controlcard,'REGULAR').gt.0.0D0) then
173 call reada(controlcard,'WEIDIS',weidis,0.1D0)
177 if (index(controlcard,'CHECKGRAD').gt.0) then
179 if (index(controlcard,'CART').gt.0) then
181 elseif (index(controlcard,'CARINT').gt.0) then
186 call reada(controlcard,'DELTA',aincr,1.0d-7)
187 elseif (index(controlcard,'THREAD').gt.0) then
189 call readi(controlcard,'THREAD',nthread,0)
190 if (nthread.gt.0) then
191 call reada(controlcard,'WEIDIS',weidis,0.1D0)
194 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
195 stop 'Error termination in Read_Control.'
197 else if (index(controlcard,'MCMA').gt.0) then
199 else if (index(controlcard,'MCEE').gt.0) then
201 else if (index(controlcard,'MULTCONF').gt.0) then
203 else if (index(controlcard,'MAP').gt.0) then
205 call readi(controlcard,'MAP',nmap,0)
206 else if (index(controlcard,'CSA').gt.0) then
208 crc else if (index(controlcard,'ZSCORE').gt.0) then
210 crc ZSCORE is rm from UNRES, modecalc=9 is available
213 cfcm else if (index(controlcard,'MCMF').gt.0) then
215 else if (index(controlcard,'SOFTREG').gt.0) then
217 else if (index(controlcard,'CHECK_BOND').gt.0) then
219 else if (index(controlcard,'TEST').gt.0) then
221 else if (index(controlcard,'MD').gt.0) then
223 else if (index(controlcard,'RE ').gt.0) then
227 lmuca=index(controlcard,'MUCA').gt.0
228 call readi(controlcard,'MUCADYN',mucadyn,0)
229 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
230 if (lmuca .and. (me.eq.king .or. .not.out1file ))
232 write (iout,*) 'MUCADYN=',mucadyn
233 write (iout,*) 'MUCASMOOTH=',muca_smooth
236 iscode=index(controlcard,'ONE_LETTER')
237 indphi=index(controlcard,'PHI')
238 indback=index(controlcard,'BACK')
239 iranconf=index(controlcard,'RAND_CONF')
240 i2ndstr=index(controlcard,'USE_SEC_PRED')
241 gradout=index(controlcard,'GRADOUT').gt.0
242 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
243 C DISTCHAINMAX become obsolete for periodic boundry condition
244 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
245 C Reading the dimensions of box in x,y,z coordinates
246 call reada(controlcard,'BOXX',boxxsize,100.0d0)
247 call reada(controlcard,'BOXY',boxysize,100.0d0)
248 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
249 c Cutoff range for interactions
250 call reada(controlcard,"R_CUT",r_cut,15.0d0)
251 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
252 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
253 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
254 if (lipthick.gt.0.0d0) then
255 bordliptop=(boxzsize+lipthick)/2.0
256 bordlipbot=bordliptop-lipthick
258 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
259 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
260 buflipbot=bordlipbot+lipbufthick
261 bufliptop=bordliptop-lipbufthick
262 if ((lipbufthick*2.0d0).gt.lipthick)
263 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
265 write(iout,*) "bordliptop=",bordliptop
266 write(iout,*) "bordlipbot=",bordlipbot
267 write(iout,*) "bufliptop=",bufliptop
268 write(iout,*) "buflipbot=",buflipbot
269 write (iout,*) "SHIELD MODE",shield_mode
270 if (TUBElog.gt.0) then
271 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
272 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
273 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
274 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
275 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
276 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
277 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
278 buftubebot=bordtubebot+tubebufthick
279 buftubetop=bordtubetop-tubebufthick
281 if (shield_mode.gt.0) then
283 C VSolvSphere the volume of solving sphere
285 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
286 C there will be no distinction between proline peptide group and normal peptide
287 C group in case of shielding parameters
288 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
289 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
290 write (iout,*) VSolvSphere,VSolvSphere_div
291 C long axis of side chain
293 long_r_sidechain(i)=vbldsc0(1,i)
294 short_r_sidechain(i)=sigma0(i)
298 if (me.eq.king .or. .not.out1file )
299 & write (iout,*) "DISTCHAINMAX",distchainmax
301 if(me.eq.king.or..not.out1file)
302 & write (iout,'(2a)') diagmeth(kdiag),
303 & ' routine used to diagonalize matrices.'
306 c--------------------------------------------------------------------------
307 subroutine read_REMDpar
311 implicit real*8 (a-h,o-z)
313 include 'COMMON.IOUNITS'
314 include 'COMMON.TIME1'
317 include 'COMMON.LANGEVIN'
319 include 'COMMON.LANGEVIN.lang0'
321 include 'COMMON.INTERACT'
322 include 'COMMON.NAMES'
324 include 'COMMON.REMD'
325 include 'COMMON.CONTROL'
326 include 'COMMON.SETUP'
328 character*320 controlcard
329 character*3200 controlcard1
330 integer iremd_m_total
332 if(me.eq.king.or..not.out1file)
333 & write (iout,*) "REMD setup"
335 call card_concat(controlcard)
336 call readi(controlcard,"NREP",nrep,3)
337 call readi(controlcard,"NSTEX",nstex,1000)
338 call reada(controlcard,"RETMIN",retmin,10.0d0)
339 call reada(controlcard,"RETMAX",retmax,1000.0d0)
340 mremdsync=(index(controlcard,'SYNC').gt.0)
341 call readi(controlcard,"NSYN",i_sync_step,100)
342 restart1file=(index(controlcard,'REST1FILE').gt.0)
343 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
344 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
345 if(max_cache_traj_use.gt.max_cache_traj)
346 & max_cache_traj_use=max_cache_traj
347 if(me.eq.king.or..not.out1file) then
348 cd if (traj1file) then
349 crc caching is in testing - NTWX is not ignored
350 cd write (iout,*) "NTWX value is ignored"
351 cd write (iout,*) " trajectory is stored to one file by master"
352 cd write (iout,*) " before exchange at NSTEX intervals"
354 write (iout,*) "NREP= ",nrep
355 write (iout,*) "NSTEX= ",nstex
356 write (iout,*) "SYNC= ",mremdsync
357 write (iout,*) "NSYN= ",i_sync_step
358 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
361 if (index(controlcard,'TLIST').gt.0) then
363 call card_concat(controlcard1)
364 read(controlcard1,*) (remd_t(i),i=1,nrep)
365 if(me.eq.king.or..not.out1file)
366 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
369 if (index(controlcard,'MLIST').gt.0) then
371 call card_concat(controlcard1)
372 read(controlcard1,*) (remd_m(i),i=1,nrep)
373 if(me.eq.king.or..not.out1file) then
374 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
377 iremd_m_total=iremd_m_total+remd_m(i)
379 write (iout,*) 'Total number of replicas ',iremd_m_total
382 if(me.eq.king.or..not.out1file)
383 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
386 c--------------------------------------------------------------------------
387 subroutine read_MDpar
391 implicit real*8 (a-h,o-z)
393 include 'COMMON.IOUNITS'
394 include 'COMMON.TIME1'
397 include 'COMMON.LANGEVIN'
399 include 'COMMON.LANGEVIN.lang0'
401 include 'COMMON.INTERACT'
402 include 'COMMON.NAMES'
404 include 'COMMON.SETUP'
405 include 'COMMON.CONTROL'
406 include 'COMMON.SPLITELE'
408 character*320 controlcard
410 call card_concat(controlcard)
411 call readi(controlcard,"NSTEP",n_timestep,1000000)
412 call readi(controlcard,"NTWE",ntwe,100)
413 call readi(controlcard,"NTWX",ntwx,1000)
414 call reada(controlcard,"DT",d_time,1.0d-1)
415 call reada(controlcard,"DVMAX",dvmax,2.0d1)
416 call reada(controlcard,"DAMAX",damax,1.0d1)
417 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
418 call readi(controlcard,"LANG",lang,0)
419 RESPA = index(controlcard,"RESPA") .gt. 0
420 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
421 ntime_split0=ntime_split
422 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
423 ntime_split0=ntime_split
424 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
425 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
426 rest = index(controlcard,"REST").gt.0
427 tbf = index(controlcard,"TBF").gt.0
428 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
429 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
430 tnh = index(controlcard,"NOSEHOOVER96").gt.0
431 if (RESPA.and.tnh)then
432 xiresp = index(controlcard,"XIRESP").gt.0
434 write(iout,*) "tnh",tnh
435 call reada(controlcard,"Q_NP",Q_np,0.1d0)
436 usampl = index(controlcard,"USAMPL").gt.0
438 mdpdb = index(controlcard,"MDPDB").gt.0
439 call reada(controlcard,"T_BATH",t_bath,300.0d0)
440 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
441 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
442 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
443 if (count_reset_moment.eq.0) count_reset_moment=1000000000
444 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
445 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
446 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
447 if (count_reset_vel.eq.0) count_reset_vel=1000000000
448 large = index(controlcard,"LARGE").gt.0
449 print_compon = index(controlcard,"PRINT_COMPON").gt.0
450 rattle = index(controlcard,"RATTLE").gt.0
451 c if performing umbrella sampling, fragments constrained are read from the fragment file
457 if(me.eq.king.or..not.out1file) then
459 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
461 write (iout,'(a)') "The units are:"
462 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
463 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
464 & " acceleration: angstrom/(48.9 fs)**2"
465 write (iout,'(a)') "energy: kcal/mol, temperature: K"
467 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
468 write (iout,'(a60,f10.5,a)')
469 & "Initial time step of numerical integration:",d_time,
471 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
473 write (iout,'(2a,i4,a)')
474 & "A-MTS algorithm used; initial time step for fast-varying",
475 & " short-range forces split into",ntime_split," steps."
476 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
477 & r_cut," lambda",rlamb
479 write (iout,'(2a,f10.5)')
480 & "Maximum acceleration threshold to reduce the time step",
481 & "/increase split number:",damax
482 write (iout,'(2a,f10.5)')
483 & "Maximum predicted energy drift to reduce the timestep",
484 & "/increase split number:",edriftmax
485 write (iout,'(a60,f10.5)')
486 & "Maximum velocity threshold to reduce velocities:",dvmax
487 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
488 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
489 if (rattle) write (iout,'(a60)')
490 & "Rattle algorithm used to constrain the virtual bonds"
494 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
495 call reada(controlcard,"RWAT",rwat,1.4d0)
496 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
497 surfarea=index(controlcard,"SURFAREA").gt.0
498 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
499 if(me.eq.king.or..not.out1file)then
500 write (iout,'(/a,$)') "Langevin dynamics calculation"
503 & " with direct integration of Langevin equations"
504 else if (lang.eq.2) then
505 write (iout,'(a/)') " with TINKER stochasic MD integrator"
506 else if (lang.eq.3) then
507 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
508 else if (lang.eq.4) then
509 write (iout,'(a/)') " in overdamped mode"
511 write (iout,'(//a,i5)')
512 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
515 write (iout,'(a60,f10.5)') "Temperature:",t_bath
516 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
517 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
518 write (iout,'(a60,f10.5)')
519 & "Scaling factor of the friction forces:",scal_fric
520 if (surfarea) write (iout,'(2a,i10,a)')
521 & "Friction coefficients will be scaled by solvent-accessible",
522 & " surface area every",reset_fricmat," steps."
524 c Calculate friction coefficients and bounds of stochastic forces
525 eta=6*pi*cPoise*etawat
526 if(me.eq.king.or..not.out1file)
527 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
529 gamp=scal_fric*(pstok+rwat)*eta
530 stdfp=dsqrt(2*Rb*t_bath/d_time)
532 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
533 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
535 if(me.eq.king.or..not.out1file)then
536 write (iout,'(/2a/)')
537 & "Radii of site types and friction coefficients and std's of",
538 & " stochastic forces of fully exposed sites"
539 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
541 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
542 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
546 if(me.eq.king.or..not.out1file)then
547 write (iout,'(a)') "Berendsen bath calculation"
548 write (iout,'(a60,f10.5)') "Temperature:",t_bath
549 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
551 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
552 & count_reset_moment," steps"
554 & write (iout,'(a,i10,a)')
555 & "Velocities will be reset at random every",count_reset_vel,
558 else if (tnp .or. tnp1 .or. tnh) then
559 if (tnp .or. tnp1) then
560 write (iout,'(a)') "Nose-Poincare bath calculation"
561 if (tnp) write (iout,'(a)')
562 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
563 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
565 write (iout,'(a)') "Nose-Hoover bath calculation"
566 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
576 WDTI(i) = 1.0*d_time/nresn
583 write (iout,'(a)') "NVT-XI-RESPA algorithm"
585 write (iout,'(a)') "NVT-XO-RESPA algorithm"
588 WDTIi(i) = 1.0*d_time/nresn/ntime_split
596 if(me.eq.king.or..not.out1file)
597 & write (iout,'(a31)') "Microcanonical mode calculation"
599 if(me.eq.king.or..not.out1file)then
600 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
602 write(iout,*) "MD running with constraints."
603 write(iout,*) "Equilibration time ", eq_time, " mtus."
604 write(iout,*) "Constraining ", nfrag," fragments."
605 write(iout,*) "Length of each fragment, weight and q0:"
607 write (iout,*) "Set of restraints #",iset
609 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
610 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
612 write(iout,*) "constraints between ", npair, "fragments."
613 write(iout,*) "constraint pairs, weights and q0:"
615 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
616 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
618 write(iout,*) "angle constraints within ", nfrag_back,
619 & "backbone fragments."
620 write(iout,*) "fragment, weights:"
622 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
623 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
624 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
627 iset=mod(kolor,nset)+1
630 if(me.eq.king.or..not.out1file)
631 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
634 c------------------------------------------------------------------------------
637 C Read molecular data.
639 implicit real*8 (a-h,o-z)
645 include 'COMMON.IOUNITS'
648 include 'COMMON.INTERACT'
649 include 'COMMON.LOCAL'
650 include 'COMMON.NAMES'
651 include 'COMMON.CHAIN'
652 include 'COMMON.FFIELD'
653 include 'COMMON.SBRIDGE'
654 include 'COMMON.HEADER'
655 include 'COMMON.CONTROL'
656 include 'COMMON.DBASE'
657 include 'COMMON.THREAD'
658 include 'COMMON.CONTACTS'
659 include 'COMMON.TORCNSTR'
660 include 'COMMON.TIME1'
661 include 'COMMON.BOUNDS'
663 include 'COMMON.SETUP'
664 include 'COMMON.SHIELD'
665 character*4 sequence(maxres)
667 double precision x(maxvar)
668 character*256 pdbfile
669 character*400 weightcard
670 character*80 weightcard_t,ucase
671 dimension itype_pdb(maxres)
672 common /pizda/ itype_pdb
673 logical seq_comp,fail
674 double precision energia(0:n_ene)
680 C Read weights of the subsequent energy terms.
681 call card_concat(weightcard)
682 call reada(weightcard,'WLONG',wlong,1.0D0)
683 call reada(weightcard,'WSC',wsc,wlong)
684 call reada(weightcard,'WSCP',wscp,wlong)
685 call reada(weightcard,'WELEC',welec,1.0D0)
686 call reada(weightcard,'WVDWPP',wvdwpp,welec)
687 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
688 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
689 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
690 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
691 call reada(weightcard,'WTURN3',wturn3,1.0D0)
692 call reada(weightcard,'WTURN4',wturn4,1.0D0)
693 call reada(weightcard,'WTURN6',wturn6,1.0D0)
694 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
695 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
696 call reada(weightcard,'WBOND',wbond,1.0D0)
697 call reada(weightcard,'WTOR',wtor,1.0D0)
698 call reada(weightcard,'WTORD',wtor_d,1.0D0)
699 call reada(weightcard,'WANG',wang,1.0D0)
700 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
701 call reada(weightcard,'SCAL14',scal14,0.4D0)
702 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
703 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
704 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
705 call reada(weightcard,'TEMP0',temp0,300.0d0)
706 call reada(weightcard,'WSHIELD',wshield,1.0d0)
707 call reada(weightcard,'WLT',wliptran,0.0D0)
708 call reada(weightcard,'WTUBE',wtube,1.0D0)
709 if (index(weightcard,'SOFT').gt.0) ipot=6
710 C 12/1/95 Added weight for the multi-body term WCORR
711 call reada(weightcard,'WCORRH',wcorr,1.0D0)
712 if (wcorr4.gt.0.0d0) wcorr=wcorr4
732 if(me.eq.king.or..not.out1file)
733 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
734 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
736 10 format (/'Energy-term weights (unscaled):'//
737 & 'WSCC= ',f10.6,' (SC-SC)'/
738 & 'WSCP= ',f10.6,' (SC-p)'/
739 & 'WELEC= ',f10.6,' (p-p electr)'/
740 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
741 & 'WBOND= ',f10.6,' (stretching)'/
742 & 'WANG= ',f10.6,' (bending)'/
743 & 'WSCLOC= ',f10.6,' (SC local)'/
744 & 'WTOR= ',f10.6,' (torsional)'/
745 & 'WTORD= ',f10.6,' (double torsional)'/
746 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
747 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
748 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
749 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
750 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
751 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
752 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
753 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
754 & 'WTURN6= ',f10.6,' (turns, 6th order)')
755 if(me.eq.king.or..not.out1file)then
756 if (wcorr4.gt.0.0d0) then
757 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
758 & 'between contact pairs of peptide groups'
759 write (iout,'(2(a,f5.3/))')
760 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
761 & 'Range of quenching the correlation terms:',2*delt_corr
762 else if (wcorr.gt.0.0d0) then
763 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
764 & 'between contact pairs of peptide groups'
766 write (iout,'(a,f8.3)')
767 & 'Scaling factor of 1,4 SC-p interactions:',scal14
768 write (iout,'(a,f8.3)')
769 & 'General scaling factor of SC-p interactions:',scalscp
771 r0_corr=cutoff_corr-delt_corr
773 aad(i,1)=scalscp*aad(i,1)
774 aad(i,2)=scalscp*aad(i,2)
775 bad(i,1)=scalscp*bad(i,1)
776 bad(i,2)=scalscp*bad(i,2)
778 call rescale_weights(t_bath)
779 if(me.eq.king.or..not.out1file)
780 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
781 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
783 22 format (/'Energy-term weights (scaled):'//
784 & 'WSCC= ',f10.6,' (SC-SC)'/
785 & 'WSCP= ',f10.6,' (SC-p)'/
786 & 'WELEC= ',f10.6,' (p-p electr)'/
787 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
788 & 'WBOND= ',f10.6,' (stretching)'/
789 & 'WANG= ',f10.6,' (bending)'/
790 & 'WSCLOC= ',f10.6,' (SC local)'/
791 & 'WTOR= ',f10.6,' (torsional)'/
792 & 'WTORD= ',f10.6,' (double torsional)'/
793 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
794 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
795 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
796 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
797 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
798 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
799 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
800 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
801 & 'WTURN6= ',f10.6,' (turns, 6th order)')
802 if(me.eq.king.or..not.out1file)
803 & write (iout,*) "Reference temperature for weights calculation:",
805 call reada(weightcard,"D0CM",d0cm,3.78d0)
806 call reada(weightcard,"AKCM",akcm,15.1d0)
807 call reada(weightcard,"AKTH",akth,11.0d0)
808 call reada(weightcard,"AKCT",akct,12.0d0)
809 call reada(weightcard,"V1SS",v1ss,-1.08d0)
810 call reada(weightcard,"V2SS",v2ss,7.61d0)
811 call reada(weightcard,"V3SS",v3ss,13.7d0)
812 call reada(weightcard,"EBR",ebr,-5.50D0)
813 call reada(weightcard,"ATRISS",atriss,0.301D0)
814 call reada(weightcard,"BTRISS",btriss,0.021D0)
815 call reada(weightcard,"CTRISS",ctriss,1.001D0)
816 call reada(weightcard,"DTRISS",dtriss,1.001D0)
817 call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
818 write (iout,*) "ATRISS=", atriss
819 write (iout,*) "BTRISS=", btriss
820 write (iout,*) "CTRISS=", ctriss
821 write (iout,*) "DTRISS=", dtriss
822 if (shield_mode.gt.0) then
824 write (iout,*) "liscale not used in shield mode"
826 write (iout,*) "lipid scaling factor", lipscale
827 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
829 dyn_ss_mask(i)=.false.
833 dyn_ssbond_ij(i,j)=1.0d300
836 call reada(weightcard,"HT",Ht,0.0D0)
838 ss_depth=ebr/wsc-0.25*eps(1,1)
839 Ht=Ht/wsc-0.25*eps(1,1)
847 if (wstrain.ne.0.0) then
848 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
853 write (iout,*) "wshield,", wshield
854 if(me.eq.king.or..not.out1file) then
855 write (iout,*) "Parameters of the SS-bond potential:"
856 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
858 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
859 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
860 write (iout,*)" HT",Ht
861 print *,'indpdb=',indpdb,' pdbref=',pdbref
863 if (indpdb.gt.0 .or. pdbref) then
864 read(inp,'(a)') pdbfile
865 if(me.eq.king.or..not.out1file)
866 & write (iout,'(2a)') 'PDB data will be read from file ',
867 & pdbfile(:ilen(pdbfile))
868 open(ipdbin,file=pdbfile,status='old',err=33)
870 33 write (iout,'(a)') 'Error opening PDB file.'
873 c write (iout,*) 'Begin reading pdb data'
876 c write (iout,*) 'Finished reading pdb data'
878 if(me.eq.king.or..not.out1file)
879 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
880 & ' nstart_sup=',nstart_sup
882 itype_pdb(i)=itype(i)
886 nct=nstart_sup+nsup-1
887 call contact(.false.,ncont_ref,icont_ref,co)
890 if(me.eq.king.or..not.out1file)
891 & write(iout,*)'Adding sidechains'
895 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
898 do while (fail.and.nsi.le.maxsi)
899 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
902 if(fail) write(iout,*)'Adding sidechain failed for res ',
903 & i,' after ',nsi,' trials'
908 if (indpdb.eq.0) then
909 C Read sequence if not taken from the pdb file.
911 c print *,'nres=',nres
912 if (iscode.gt.0) then
913 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
915 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
917 C Convert sequence to numeric code
919 itype(i)=rescode(i,sequence(i),iscode)
921 C Assign initial virtual bond lengths
927 vbld(i+nres)=dsc(iabs(itype(i)))
928 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
929 c write (iout,*) "i",i," itype",itype(i),
930 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
934 c print '(20i4)',(itype(i),i=1,nres)
937 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
939 if (itype(i).eq.ntyp1) then
943 else if (iabs(itype(i+1)).ne.20) then
945 else if (iabs(itype(i)).ne.20) then
952 if(me.eq.king.or..not.out1file)then
953 write (iout,*) "ITEL"
955 write (iout,*) i,itype(i),itel(i)
957 print *,'Call Read_Bridge.'
960 C 8/13/98 Set limits to generating the dihedral angles
965 read (inp,*) ndih_constr
966 if (ndih_constr.gt.0) then
968 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
970 if(me.eq.king.or..not.out1file)then
972 & 'There are',ndih_constr,' constraints on phi angles.'
974 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
979 phi0(i)=deg2rad*phi0(i)
980 drange(i)=deg2rad*drange(i)
982 C if(me.eq.king.or..not.out1file)
983 C & write (iout,*) 'FTORS',ftors
986 phibound(1,ii) = phi0(i)-drange(i)
987 phibound(2,ii) = phi0(i)+drange(i)
990 C first setting the theta boundaries to 0 to pi
991 C this mean that there is no energy penalty for any angle occuring this can be applied
992 C for generate random conformation but is not implemented in this way
997 C begin reading theta constrains this is quartic constrains allowing to
998 C have smooth second derivative
999 if (with_theta_constr) then
1000 C with_theta_constr is keyword allowing for occurance of theta constrains
1001 read (inp,*) ntheta_constr
1002 C ntheta_constr is the number of theta constrains
1003 if (ntheta_constr.gt.0) then
1004 C read (inp,*) ftors
1005 read (inp,*) (itheta_constr(i),theta_constr0(i),
1006 & theta_drange(i),for_thet_constr(i),
1007 & i=1,ntheta_constr)
1008 C the above code reads from 1 to ntheta_constr
1009 C itheta_constr(i) residue i for which is theta_constr
1010 C theta_constr0 the global minimum value
1011 C theta_drange is range for which there is no energy penalty
1012 C for_thet_constr is the force constant for quartic energy penalty
1014 if(me.eq.king.or..not.out1file)then
1016 & 'There are',ntheta_constr,' constraints on phi angles.'
1017 do i=1,ntheta_constr
1018 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1020 & for_thet_constr(i)
1023 do i=1,ntheta_constr
1024 theta_constr0(i)=deg2rad*theta_constr0(i)
1025 theta_drange(i)=deg2rad*theta_drange(i)
1027 C if(me.eq.king.or..not.out1file)
1028 C & write (iout,*) 'FTORS',ftors
1029 C do i=1,ntheta_constr
1030 C ii = itheta_constr(i)
1031 C thetabound(1,ii) = phi0(i)-drange(i)
1032 C thetabound(2,ii) = phi0(i)+drange(i)
1034 endif ! ntheta_constr.gt.0
1035 endif! with_theta_constr
1037 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1038 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1041 if (me.eq.king) then
1043 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1045 write (iout,'(a3,i5,2f10.1)')
1046 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1052 cd print *,'NNT=',NNT,' NCT=',NCT
1053 if (itype(1).eq.ntyp1) nnt=2
1054 if (itype(nres).eq.ntyp1) nct=nct-1
1056 if(me.eq.king.or..not.out1file)
1057 & write (iout,'(a,i3)') 'nsup=',nsup
1059 if (nsup.le.(nct-nnt+1)) then
1060 do i=0,nct-nnt+1-nsup
1061 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1067 & 'Error - sequences to be superposed do not match.'
1070 do i=0,nsup-(nct-nnt+1)
1071 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1073 nstart_sup=nstart_sup+i
1079 & 'Error - sequences to be superposed do not match.'
1082 if (nsup.eq.0) nsup=nct-nnt
1083 if (nstart_sup.eq.0) nstart_sup=nnt
1084 if (nstart_seq.eq.0) nstart_seq=nnt
1085 if(me.eq.king.or..not.out1file)
1086 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1087 & ' nstart_seq=',nstart_seq
1089 c--- Zscore rms -------
1090 if (nz_start.eq.0) nz_start=nnt
1091 if (nz_end.eq.0 .and. nsup.gt.0) then
1093 else if (nz_end.eq.0) then
1096 if(me.eq.king.or..not.out1file)then
1097 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1098 write (iout,*) 'IZ_SC=',iz_sc
1100 c----------------------
1103 if (.not.pdbref) then
1104 call read_angles(inp,*38)
1106 38 write (iout,'(a)') 'Error reading reference structure.'
1108 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1109 stop 'Error reading reference structure'
1111 39 call chainbuild_extconf
1113 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1120 cref(j,i,kkk)=c(j,i)
1123 call contact(.true.,ncont_ref,icont_ref,co)
1127 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1129 if (constr_dist.gt.0) call read_dist_constr
1130 write (iout,*) "After read_dist_constr nhpb",nhpb
1131 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1133 if(me.eq.king.or..not.out1file)
1134 & write (iout,*) 'Contact order:',co
1136 if(me.eq.king.or..not.out1file)
1137 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1140 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1142 if(me.eq.king.or..not.out1file)
1143 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1144 & icont_ref(1,i),' ',
1145 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1149 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1150 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1151 & modecalc.ne.10) then
1152 C If input structure hasn't been supplied from the PDB file read or generate
1154 if (iranconf.eq.0 .and. .not. extconf) then
1155 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1156 & write (iout,'(a)') 'Initial geometry will be read in.'
1158 read(inp,'(8f10.5)',end=36,err=36)
1159 & ((c(l,k),l=1,3),k=1,nres),
1160 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1161 write (iout,*) "Exit READ_CART"
1162 write (iout,'(8f10.5)')
1163 & ((c(l,k),l=1,3),k=1,nres),
1164 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1165 call int_from_cart1(.true.)
1166 write (iout,*) "Finish INT_TO_CART"
1169 dc(j,i)=c(j,i+1)-c(j,i)
1170 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1174 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1176 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1177 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1183 call read_angles(inp,*36)
1184 call chainbuild_extconf
1187 36 write (iout,'(a)') 'Error reading angle file.'
1189 call mpi_finalize( MPI_COMM_WORLD,IERR )
1191 stop 'Error reading angle file.'
1193 else if (extconf) then
1194 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1195 & write (iout,'(a)') 'Extended chain initial geometry.'
1197 theta(i)=90d0*deg2rad
1200 phi(i)=180d0*deg2rad
1203 alph(i)=110d0*deg2rad
1206 omeg(i)=-120d0*deg2rad
1207 if (itype(i).le.0) omeg(i)=-omeg(i)
1209 call chainbuild_extconf
1211 if(me.eq.king.or..not.out1file)
1212 & write (iout,'(a)') 'Random-generated initial geometry.'
1216 if (me.eq.king .or. fg_rank.eq.0 .and. (
1217 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1221 call gen_rand_conf(itmp,*30)
1223 30 write (iout,*) 'Failed to generate random conformation',
1224 & ', itrial=',itrial
1225 write (*,*) 'Processor:',me,
1226 & ' Failed to generate random conformation',
1236 write (iout,'(a,i3,a)') 'Processor:',me,
1237 & ' error in generating random conformation.'
1238 write (*,'(a,i3,a)') 'Processor:',me,
1239 & ' error in generating random conformation.'
1242 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1247 & ' error in generating random conformation.'
1252 elseif (modecalc.eq.4) then
1253 read (inp,'(a)') intinname
1254 open (intin,file=intinname,status='old',err=333)
1255 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1256 & write (iout,'(a)') 'intinname',intinname
1257 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1259 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1261 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1263 stop 'Error opening angle file.'
1267 C Generate distance constraints, if the PDB structure is to be regularized.
1268 if (nthread.gt.0) then
1269 call read_threadbase
1272 if (me.eq.king .or. .not. out1file)
1274 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1275 write (iout,'(/a,i3,a)')
1276 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1277 write (iout,'(20i4)') (iss(i),i=1,ns)
1279 write(iout,*)"Running with dynamic disulfide-bond formation"
1281 write (iout,'(/a/)') 'Pre-formed links are:'
1287 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1288 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1294 if (ns.gt.0.and.dyn_ss) then
1298 forcon(i-nss)=forcon(i)
1305 dyn_ss_mask(iss(i))=.true.
1308 if (i2ndstr.gt.0) call secstrp2dihc
1309 c call geom_to_var(nvar,x)
1310 c call etotal(energia(0))
1311 c call enerprint(energia(0))
1312 c call briefout(0,etot)
1314 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1315 cd write (iout,'(a)') 'Variable list:'
1316 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1318 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1319 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1320 & 'Processor',myrank,': end reading molecular data.'
1325 c--------------------------------------------------------------------------
1326 logical function seq_comp(itypea,itypeb,length)
1328 integer length,itypea(length),itypeb(length)
1331 if (itypea(i).ne.itypeb(i)) then
1339 c-----------------------------------------------------------------------------
1340 subroutine read_bridge
1341 C Read information about disulfide bridges.
1342 implicit real*8 (a-h,o-z)
1343 include 'DIMENSIONS'
1347 include 'COMMON.IOUNITS'
1348 include 'COMMON.GEO'
1349 include 'COMMON.VAR'
1350 include 'COMMON.INTERACT'
1351 include 'COMMON.LOCAL'
1352 include 'COMMON.NAMES'
1353 include 'COMMON.CHAIN'
1354 include 'COMMON.FFIELD'
1355 include 'COMMON.SBRIDGE'
1356 include 'COMMON.HEADER'
1357 include 'COMMON.CONTROL'
1358 include 'COMMON.DBASE'
1359 include 'COMMON.THREAD'
1360 include 'COMMON.TIME1'
1361 include 'COMMON.SETUP'
1362 C Read bridging residues.
1363 read (inp,*) ns,(iss(i),i=1,ns)
1365 if(me.eq.king.or..not.out1file)
1366 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1367 C Check whether the specified bridging residues are cystines.
1369 if (itype(iss(i)).ne.1) then
1370 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1371 & 'Do you REALLY think that the residue ',
1372 & restyp(itype(iss(i))),i,
1373 & ' can form a disulfide bridge?!!!'
1374 write (*,'(2a,i3,a)')
1375 & 'Do you REALLY think that the residue ',
1376 & restyp(itype(iss(i))),i,
1377 & ' can form a disulfide bridge?!!!'
1379 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1384 C Read preformed bridges.
1386 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1388 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1391 C Check if the residues involved in bridges are in the specified list of
1392 C bridging residues.
1395 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1396 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1397 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1398 & ' contains residues present in other pairs.'
1399 write (*,'(a,i3,a)') 'Disulfide pair',i,
1400 & ' contains residues present in other pairs.'
1402 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1408 if (ihpb(i).eq.iss(j)) goto 10
1410 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1413 if (jhpb(i).eq.iss(j)) goto 20
1415 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1421 ihpb(i)=ihpb(i)+nres
1422 jhpb(i)=jhpb(i)+nres
1428 c----------------------------------------------------------------------------
1429 subroutine read_x(kanal,*)
1430 implicit real*8 (a-h,o-z)
1431 include 'DIMENSIONS'
1432 include 'COMMON.GEO'
1433 include 'COMMON.VAR'
1434 include 'COMMON.CHAIN'
1435 include 'COMMON.IOUNITS'
1436 include 'COMMON.CONTROL'
1437 include 'COMMON.LOCAL'
1438 include 'COMMON.INTERACT'
1439 c Read coordinates from input
1441 read(kanal,'(8f10.5)',end=10,err=10)
1442 & ((c(l,k),l=1,3),k=1,nres),
1443 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1446 c(j,2*nres)=c(j,nres)
1448 call int_from_cart1(.false.)
1451 dc(j,i)=c(j,i+1)-c(j,i)
1452 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1456 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1458 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1459 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1467 c----------------------------------------------------------------------------
1468 subroutine read_threadbase
1469 implicit real*8 (a-h,o-z)
1470 include 'DIMENSIONS'
1471 include 'COMMON.IOUNITS'
1472 include 'COMMON.GEO'
1473 include 'COMMON.VAR'
1474 include 'COMMON.INTERACT'
1475 include 'COMMON.LOCAL'
1476 include 'COMMON.NAMES'
1477 include 'COMMON.CHAIN'
1478 include 'COMMON.FFIELD'
1479 include 'COMMON.SBRIDGE'
1480 include 'COMMON.HEADER'
1481 include 'COMMON.CONTROL'
1482 include 'COMMON.DBASE'
1483 include 'COMMON.THREAD'
1484 include 'COMMON.TIME1'
1485 C Read pattern database for threading.
1486 read (icbase,*) nseq
1488 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1489 & nres_base(2,i),nres_base(3,i)
1490 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1492 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1493 c & nres_base(2,i),nres_base(3,i)
1494 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1498 if (weidis.eq.0.0D0) weidis=0.1D0
1507 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1508 write (iout,'(a,i5)') 'nexcl: ',nexcl
1509 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1512 c------------------------------------------------------------------------------
1513 subroutine setup_var
1514 implicit real*8 (a-h,o-z)
1515 include 'DIMENSIONS'
1516 include 'COMMON.IOUNITS'
1517 include 'COMMON.GEO'
1518 include 'COMMON.VAR'
1519 include 'COMMON.INTERACT'
1520 include 'COMMON.LOCAL'
1521 include 'COMMON.NAMES'
1522 include 'COMMON.CHAIN'
1523 include 'COMMON.FFIELD'
1524 include 'COMMON.SBRIDGE'
1525 include 'COMMON.HEADER'
1526 include 'COMMON.CONTROL'
1527 include 'COMMON.DBASE'
1528 include 'COMMON.THREAD'
1529 include 'COMMON.TIME1'
1530 C Set up variable list.
1536 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1538 ialph(i,1)=nvar+nside
1542 if (indphi.gt.0) then
1544 else if (indback.gt.0) then
1549 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1552 c----------------------------------------------------------------------------
1553 subroutine gen_dist_constr
1554 C Generate CA distance constraints.
1555 implicit real*8 (a-h,o-z)
1556 include 'DIMENSIONS'
1557 include 'COMMON.IOUNITS'
1558 include 'COMMON.GEO'
1559 include 'COMMON.VAR'
1560 include 'COMMON.INTERACT'
1561 include 'COMMON.LOCAL'
1562 include 'COMMON.NAMES'
1563 include 'COMMON.CHAIN'
1564 include 'COMMON.FFIELD'
1565 include 'COMMON.SBRIDGE'
1566 include 'COMMON.HEADER'
1567 include 'COMMON.CONTROL'
1568 include 'COMMON.DBASE'
1569 include 'COMMON.THREAD'
1570 include 'COMMON.TIME1'
1571 dimension itype_pdb(maxres)
1572 common /pizda/ itype_pdb
1574 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1575 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1576 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1578 if (constr_dist.eq.11) then
1579 do i=nstart_sup,nstart_sup+nsup-1
1580 do j=i+2,nstart_sup+nsup-1
1582 if (distance.le.15.0) then
1584 ihpb(nhpb)=i+nstart_seq-nstart_sup
1585 jhpb(nhpb)=j+nstart_seq-nstart_sup
1586 forcon(nhpb)=sqrt(0.04*distance)
1587 fordepth(nhpb)=sqrt(40.0/distance)
1588 dhpb(nhpb)=distance-0.1d0
1589 dhpb1(nhpb)=distance+0.1d0
1592 if (.not.out1file .or. me.eq.king)
1593 & write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1594 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1596 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1597 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1603 do i=nstart_sup,nstart_sup+nsup-1
1604 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1605 cd & ' seq_pdb', restyp(itype_pdb(i))
1606 do j=i+2,nstart_sup+nsup-1
1608 ihpb(nhpb)=i+nstart_seq-nstart_sup
1609 jhpb(nhpb)=j+nstart_seq-nstart_sup
1611 dhpb(nhpb)=dist(i,j)
1615 cd write (iout,'(a)') 'Distance constraints:'
1620 cd if (ii.gt.nres) then
1625 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1626 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1627 cd & dhpb(i),forcon(i)
1631 c----------------------------------------------------------------------------
1633 implicit real*8 (a-h,o-z)
1634 include 'DIMENSIONS'
1635 include 'COMMON.MAP'
1636 include 'COMMON.IOUNITS'
1637 character*3 angid(4) /'THE','PHI','ALP','OME'/
1638 character*80 mapcard,ucase
1640 read (inp,'(a)') mapcard
1641 mapcard=ucase(mapcard)
1642 if (index(mapcard,'PHI').gt.0) then
1644 else if (index(mapcard,'THE').gt.0) then
1646 else if (index(mapcard,'ALP').gt.0) then
1648 else if (index(mapcard,'OME').gt.0) then
1651 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1652 stop 'Error - illegal variable spec in MAP card.'
1654 call readi (mapcard,'RES1',res1(imap),0)
1655 call readi (mapcard,'RES2',res2(imap),0)
1656 if (res1(imap).eq.0) then
1657 res1(imap)=res2(imap)
1658 else if (res2(imap).eq.0) then
1659 res2(imap)=res1(imap)
1661 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1663 & 'Error - illegal definition of variable group in MAP.'
1664 stop 'Error - illegal definition of variable group in MAP.'
1666 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1667 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1668 call readi(mapcard,'NSTEP',nstep(imap),0)
1669 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1671 & 'Illegal boundary and/or step size specification in MAP.'
1672 stop 'Illegal boundary and/or step size specification in MAP.'
1677 c----------------------------------------------------------------------------
1679 implicit real*8 (a-h,o-z)
1680 include 'DIMENSIONS'
1681 include 'COMMON.IOUNITS'
1682 include 'COMMON.GEO'
1683 include 'COMMON.CSA'
1684 include 'COMMON.BANK'
1685 include 'COMMON.CONTROL'
1687 character*620 mcmcard
1688 call card_concat(mcmcard)
1690 call readi(mcmcard,'NCONF',nconf,50)
1691 call readi(mcmcard,'NADD',nadd,0)
1692 call readi(mcmcard,'JSTART',jstart,1)
1693 call readi(mcmcard,'JEND',jend,1)
1694 call readi(mcmcard,'NSTMAX',nstmax,500000)
1695 call readi(mcmcard,'N0',n0,1)
1696 call readi(mcmcard,'N1',n1,6)
1697 call readi(mcmcard,'N2',n2,4)
1698 call readi(mcmcard,'N3',n3,0)
1699 call readi(mcmcard,'N4',n4,0)
1700 call readi(mcmcard,'N5',n5,0)
1701 call readi(mcmcard,'N6',n6,10)
1702 call readi(mcmcard,'N7',n7,0)
1703 call readi(mcmcard,'N8',n8,0)
1704 call readi(mcmcard,'N9',n9,0)
1705 call readi(mcmcard,'N14',n14,0)
1706 call readi(mcmcard,'N15',n15,0)
1707 call readi(mcmcard,'N16',n16,0)
1708 call readi(mcmcard,'N17',n17,0)
1709 call readi(mcmcard,'N18',n18,0)
1711 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1713 call readi(mcmcard,'NDIFF',ndiff,2)
1714 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1715 call readi(mcmcard,'IS1',is1,1)
1716 call readi(mcmcard,'IS2',is2,8)
1717 call readi(mcmcard,'NRAN0',nran0,4)
1718 call readi(mcmcard,'NRAN1',nran1,2)
1719 call readi(mcmcard,'IRR',irr,1)
1720 call readi(mcmcard,'NSEED',nseed,20)
1721 call readi(mcmcard,'NTOTAL',ntotal,10000)
1722 call reada(mcmcard,'CUT1',cut1,2.0d0)
1723 call reada(mcmcard,'CUT2',cut2,5.0d0)
1724 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1725 call readi(mcmcard,'ICMAX',icmax,3)
1726 call readi(mcmcard,'IRESTART',irestart,0)
1727 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1730 call reada(mcmcard,'DELE',dele,20.0d0)
1731 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1732 call readi(mcmcard,'IREF',iref,0)
1733 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1734 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1735 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1736 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1737 write (iout,*) "NCONF_IN",nconf_in
1740 c----------------------------------------------------------------------------
1741 cfmc subroutine mcmfread
1742 cfmc implicit real*8 (a-h,o-z)
1743 cfmc include 'DIMENSIONS'
1744 cfmc include 'COMMON.MCMF'
1745 cfmc include 'COMMON.IOUNITS'
1746 cfmc include 'COMMON.GEO'
1747 cfmc character*80 ucase
1748 cfmc character*620 mcmcard
1749 cfmc call card_concat(mcmcard)
1751 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1752 cfmc write(iout,*)'MAXRANT=',maxrant
1753 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1754 cfmc write(iout,*)'MAXFAM=',maxfam
1755 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1756 cfmc write(iout,*)'NNET1=',nnet1
1757 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1758 cfmc write(iout,*)'NNET2=',nnet2
1759 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1760 cfmc write(iout,*)'NNET3=',nnet3
1761 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1762 cfmc write(iout,*)'ILASTT=',ilastt
1763 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1764 cfmc write(iout,*)'MAXSTR=',maxstr
1765 cfmc maxstr_f=maxstr/maxfam
1766 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1767 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1768 cfmc write(iout,*)'NMCMF=',nmcmf
1769 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1770 cfmc write(iout,*)'IFOCUS=',ifocus
1771 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1772 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1773 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1774 cfmc write(iout,*)'INTPRT=',intprt
1775 cfmc call readi(mcmcard,'IPRT',iprt,100)
1776 cfmc write(iout,*)'IPRT=',iprt
1777 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1778 cfmc write(iout,*)'IMAXTR=',imaxtr
1779 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1780 cfmc write(iout,*)'MAXEVEN=',maxeven
1781 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1782 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1783 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1784 cfmc write(iout,*)'INIMIN=',inimin
1785 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1786 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1787 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1788 cfmc write(iout,*)'NTHREAD=',nthread
1789 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1790 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1791 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1792 cfmc write(iout,*)'MAXPERT=',maxpert
1793 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1794 cfmc write(iout,*)'IRMSD=',irmsd
1795 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1796 cfmc write(iout,*)'DENEMIN=',denemin
1797 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1798 cfmc write(iout,*)'RCUT1S=',rcut1s
1799 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1800 cfmc write(iout,*)'RCUT1E=',rcut1e
1801 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1802 cfmc write(iout,*)'RCUT2S=',rcut2s
1803 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1804 cfmc write(iout,*)'RCUT2E=',rcut2e
1805 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1806 cfmc write(iout,*)'DPERT1=',d_pert1
1807 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1808 cfmc write(iout,*)'DPERT1A=',d_pert1a
1809 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1810 cfmc write(iout,*)'DPERT2=',d_pert2
1811 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1812 cfmc write(iout,*)'DPERT2A=',d_pert2a
1813 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1814 cfmc write(iout,*)'DPERT2B=',d_pert2b
1815 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1816 cfmc write(iout,*)'DPERT2C=',d_pert2c
1817 cfmc d_pert1=deg2rad*d_pert1
1818 cfmc d_pert1a=deg2rad*d_pert1a
1819 cfmc d_pert2=deg2rad*d_pert2
1820 cfmc d_pert2a=deg2rad*d_pert2a
1821 cfmc d_pert2b=deg2rad*d_pert2b
1822 cfmc d_pert2c=deg2rad*d_pert2c
1823 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1824 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1825 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1826 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1827 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1828 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1829 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1830 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1831 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1832 cfmc write(iout,*)'RCUTINI=',rcutini
1833 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1834 cfmc write(iout,*)'GRAT=',grat
1835 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1836 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1840 c----------------------------------------------------------------------------
1842 implicit real*8 (a-h,o-z)
1843 include 'DIMENSIONS'
1844 include 'COMMON.MCM'
1845 include 'COMMON.MCE'
1846 include 'COMMON.IOUNITS'
1848 character*320 mcmcard
1849 call card_concat(mcmcard)
1850 call readi(mcmcard,'MAXACC',maxacc,100)
1851 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1852 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1853 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1854 call readi(mcmcard,'MAXREPM',maxrepm,200)
1855 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1856 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1857 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1858 call reada(mcmcard,'E_UP',e_up,5.0D0)
1859 call reada(mcmcard,'DELTE',delte,0.1D0)
1860 call readi(mcmcard,'NSWEEP',nsweep,5)
1861 call readi(mcmcard,'NSTEPH',nsteph,0)
1862 call readi(mcmcard,'NSTEPC',nstepc,0)
1863 call reada(mcmcard,'TMIN',tmin,298.0D0)
1864 call reada(mcmcard,'TMAX',tmax,298.0D0)
1865 call readi(mcmcard,'NWINDOW',nwindow,0)
1866 call readi(mcmcard,'PRINT_MC',print_mc,0)
1867 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1868 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1869 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1870 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1871 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1872 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1873 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1874 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1875 if (nwindow.gt.0) then
1876 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1878 winlen(i)=winend(i)-winstart(i)+1
1881 if (tmax.lt.tmin) tmax=tmin
1882 if (tmax.eq.tmin) then
1886 if (nstepc.gt.0 .and. nsteph.gt.0) then
1887 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1888 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1890 C Probabilities of different move types
1891 sumpro_type(0)=0.0D0
1892 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1893 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1894 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1895 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1896 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1897 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1898 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1900 print *,'i',i,' sumprotype',sumpro_type(i)
1901 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1902 print *,'i',i,' sumprotype',sumpro_type(i)
1906 c----------------------------------------------------------------------------
1907 subroutine read_minim
1908 implicit real*8 (a-h,o-z)
1909 include 'DIMENSIONS'
1910 include 'COMMON.MINIM'
1911 include 'COMMON.IOUNITS'
1913 character*320 minimcard
1914 call card_concat(minimcard)
1915 call readi(minimcard,'MAXMIN',maxmin,2000)
1916 call readi(minimcard,'MAXFUN',maxfun,5000)
1917 call readi(minimcard,'MINMIN',minmin,maxmin)
1918 call readi(minimcard,'MINFUN',minfun,maxmin)
1919 call reada(minimcard,'TOLF',tolf,1.0D-2)
1920 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1921 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1922 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1923 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1924 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1925 & 'Options in energy minimization:'
1926 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1927 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1928 & 'MinMin:',MinMin,' MinFun:',MinFun,
1929 & ' TolF:',TolF,' RTolF:',RTolF
1932 c----------------------------------------------------------------------------
1933 subroutine read_angles(kanal,*)
1934 implicit real*8 (a-h,o-z)
1935 include 'DIMENSIONS'
1936 include 'COMMON.GEO'
1937 include 'COMMON.VAR'
1938 include 'COMMON.CHAIN'
1939 include 'COMMON.IOUNITS'
1940 include 'COMMON.CONTROL'
1941 c Read angles from input
1943 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1944 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1945 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1946 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1949 c 9/7/01 avoid 180 deg valence angle
1950 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1952 theta(i)=deg2rad*theta(i)
1953 phi(i)=deg2rad*phi(i)
1954 alph(i)=deg2rad*alph(i)
1955 omeg(i)=deg2rad*omeg(i)
1960 c----------------------------------------------------------------------------
1961 subroutine reada(rekord,lancuch,wartosc,default)
1963 character*(*) rekord,lancuch
1964 double precision wartosc,default
1967 iread=index(rekord,lancuch)
1968 if (iread.eq.0) then
1972 iread=iread+ilen(lancuch)+1
1973 read (rekord(iread:),*,err=10,end=10) wartosc
1978 c----------------------------------------------------------------------------
1979 subroutine readi(rekord,lancuch,wartosc,default)
1981 character*(*) rekord,lancuch
1982 integer wartosc,default
1985 iread=index(rekord,lancuch)
1986 if (iread.eq.0) then
1990 iread=iread+ilen(lancuch)+1
1991 read (rekord(iread:),*,err=10,end=10) wartosc
1996 c----------------------------------------------------------------------------
1997 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2000 integer tablica(dim),default
2001 character*(*) rekord,lancuch
2008 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2009 if (iread.eq.0) return
2010 iread=iread+ilen(lancuch)+1
2012 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2015 c----------------------------------------------------------------------------
2016 subroutine multreada(rekord,lancuch,tablica,dim,default)
2019 double precision tablica(dim),default
2020 character*(*) rekord,lancuch
2027 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2028 if (iread.eq.0) return
2029 iread=iread+ilen(lancuch)+1
2030 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2033 c----------------------------------------------------------------------------
2034 subroutine openunits
2035 implicit real*8 (a-h,o-z)
2036 include 'DIMENSIONS'
2039 character*16 form,nodename
2042 include 'COMMON.SETUP'
2043 include 'COMMON.IOUNITS'
2045 include 'COMMON.CONTROL'
2046 integer lenpre,lenpot,ilen,lentmp
2048 character*3 out1file_text,ucase
2051 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2052 call getenv_loc("PREFIX",prefix)
2054 call getenv_loc("POT",pot)
2055 call getenv_loc("DIRTMP",tmpdir)
2056 call getenv_loc("CURDIR",curdir)
2057 call getenv_loc("OUT1FILE",out1file_text)
2058 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2059 out1file_text=ucase(out1file_text)
2060 if (out1file_text(1:1).eq."Y") then
2063 out1file=fg_rank.gt.0
2068 if (lentmp.gt.0) then
2069 write (*,'(80(1h!))')
2070 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2071 write (*,'(80(1h!))')
2072 write (*,*)"All output files will be on node /tmp directory."
2074 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2075 if (me.eq.king) then
2076 write (*,*) "The master node is ",nodename
2077 else if (fg_rank.eq.0) then
2078 write (*,*) "I am the CG slave node ",nodename
2080 write (*,*) "I am the FG slave node ",nodename
2083 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2084 lenpre = lentmp+lenpre+1
2086 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2087 C Get the names and open the input files
2088 #if defined(WINIFL) || defined(WINPGI)
2089 open(1,file=pref_orig(:ilen(pref_orig))//
2090 & '.inp',status='old',readonly,shared)
2091 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2092 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2093 C Get parameter filenames and open the parameter files.
2094 call getenv_loc('BONDPAR',bondname)
2095 open (ibond,file=bondname,status='old',readonly,shared)
2096 call getenv_loc('THETPAR',thetname)
2097 open (ithep,file=thetname,status='old',readonly,shared)
2098 call getenv_loc('ROTPAR',rotname)
2099 open (irotam,file=rotname,status='old',readonly,shared)
2100 call getenv_loc('TORPAR',torname)
2101 open (itorp,file=torname,status='old',readonly,shared)
2102 call getenv_loc('TORDPAR',tordname)
2103 open (itordp,file=tordname,status='old',readonly,shared)
2104 call getenv_loc('TORKCC',torkccname)
2105 open (itorkcc,file=torkccname,status='old',readonly,shared)
2106 call getenv_loc('THETKCC',thetkccname)
2107 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2108 call getenv_loc('FOURIER',fouriername)
2109 open (ifourier,file=fouriername,status='old',readonly,shared)
2110 call getenv_loc('ELEPAR',elename)
2111 open (ielep,file=elename,status='old',readonly,shared)
2112 call getenv_loc('SIDEPAR',sidename)
2113 open (isidep,file=sidename,status='old',readonly,shared)
2114 call getenv_loc('LIPTRANPAR',liptranname)
2115 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2116 #elif (defined CRAY) || (defined AIX)
2117 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2119 c print *,"Processor",myrank," opened file 1"
2120 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2121 c print *,"Processor",myrank," opened file 9"
2122 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2123 C Get parameter filenames and open the parameter files.
2124 call getenv_loc('BONDPAR',bondname)
2125 open (ibond,file=bondname,status='old',action='read')
2126 c print *,"Processor",myrank," opened file IBOND"
2127 call getenv_loc('THETPAR',thetname)
2128 open (ithep,file=thetname,status='old',action='read')
2129 c print *,"Processor",myrank," opened file ITHEP"
2130 call getenv_loc('ROTPAR',rotname)
2131 open (irotam,file=rotname,status='old',action='read')
2132 c print *,"Processor",myrank," opened file IROTAM"
2133 call getenv_loc('TORPAR',torname)
2134 open (itorp,file=torname,status='old',action='read')
2135 c print *,"Processor",myrank," opened file ITORP"
2136 call getenv_loc('TORDPAR',tordname)
2137 open (itordp,file=tordname,status='old',action='read')
2138 call getenv_loc('TORKCC',torkccname)
2139 open (itorkcc,file=torkccname,status='old',action='read')
2140 call getenv_loc('THETKCC',thetkccname)
2141 open (ithetkcc,file=thetkccname,status='old',action='read')
2142 c print *,"Processor",myrank," opened file ITORDP"
2143 call getenv_loc('SCCORPAR',sccorname)
2144 open (isccor,file=sccorname,status='old',action='read')
2145 c print *,"Processor",myrank," opened file ISCCOR"
2146 call getenv_loc('FOURIER',fouriername)
2147 open (ifourier,file=fouriername,status='old',action='read')
2148 c print *,"Processor",myrank," opened file IFOURIER"
2149 call getenv_loc('ELEPAR',elename)
2150 open (ielep,file=elename,status='old',action='read')
2151 c print *,"Processor",myrank," opened file IELEP"
2152 call getenv_loc('SIDEPAR',sidename)
2153 open (isidep,file=sidename,status='old',action='read')
2154 call getenv_loc('LIPTRANPAR',liptranname)
2155 open (iliptranpar,file=liptranname,status='old',action='read')
2156 c print *,"Processor",myrank," opened file ISIDEP"
2157 c print *,"Processor",myrank," opened parameter files"
2159 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2160 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2161 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2162 C Get parameter filenames and open the parameter files.
2163 call getenv_loc('BONDPAR',bondname)
2164 open (ibond,file=bondname,status='old')
2165 call getenv_loc('THETPAR',thetname)
2166 open (ithep,file=thetname,status='old')
2167 call getenv_loc('ROTPAR',rotname)
2168 open (irotam,file=rotname,status='old')
2169 call getenv_loc('TORPAR',torname)
2170 open (itorp,file=torname,status='old')
2171 call getenv_loc('TORDPAR',tordname)
2172 open (itordp,file=tordname,status='old')
2173 call getenv_loc('TORKCC',torkccname)
2174 open (itorkcc,file=torkccname,status='old')
2175 call getenv_loc('THETKCC',thetkccname)
2176 open (ithetkcc,file=thetkccname,status='old')
2177 call getenv_loc('SCCORPAR',sccorname)
2178 open (isccor,file=sccorname,status='old')
2179 call getenv_loc('FOURIER',fouriername)
2180 open (ifourier,file=fouriername,status='old')
2181 call getenv_loc('ELEPAR',elename)
2182 open (ielep,file=elename,status='old')
2183 call getenv_loc('SIDEPAR',sidename)
2184 open (isidep,file=sidename,status='old')
2185 call getenv_loc('LIPTRANPAR',liptranname)
2186 open (iliptranpar,file=liptranname,status='old')
2188 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2190 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2191 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2192 C Get parameter filenames and open the parameter files.
2193 call getenv_loc('BONDPAR',bondname)
2194 open (ibond,file=bondname,status='old',readonly)
2195 call getenv_loc('THETPAR',thetname)
2196 open (ithep,file=thetname,status='old',readonly)
2197 call getenv_loc('ROTPAR',rotname)
2198 open (irotam,file=rotname,status='old',readonly)
2199 call getenv_loc('TORPAR',torname)
2200 open (itorp,file=torname,status='old',readonly)
2201 call getenv_loc('TORDPAR',tordname)
2202 open (itordp,file=tordname,status='old',readonly)
2203 call getenv_loc('TORKCC',torkccname)
2204 open (itorkcc,file=torkccname,status='old',readonly)
2205 call getenv_loc('THETKCC',thetkccname)
2206 open (ithetkcc,file=thetkccname,status='old',readonly)
2207 call getenv_loc('SCCORPAR',sccorname)
2208 open (isccor,file=sccorname,status='old',readonly)
2210 call getenv_loc('THETPARPDB',thetname_pdb)
2211 print *,"thetname_pdb ",thetname_pdb
2212 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2213 print *,ithep_pdb," opened"
2215 call getenv_loc('FOURIER',fouriername)
2216 open (ifourier,file=fouriername,status='old',readonly)
2217 call getenv_loc('ELEPAR',elename)
2218 open (ielep,file=elename,status='old',readonly)
2219 call getenv_loc('SIDEPAR',sidename)
2220 open (isidep,file=sidename,status='old',readonly)
2221 call getenv_loc('LIPTRANPAR',liptranname)
2222 open (iliptranpar,file=liptranname,status='old',action='read')
2224 call getenv_loc('ROTPARPDB',rotname_pdb)
2225 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2228 call getenv_loc('TUBEPAR',tubename)
2229 open (itube,file=tubename,status='old',readonly)
2232 C 8/9/01 In the newest version SCp interaction constants are read from a file
2233 C Use -DOLDSCP to use hard-coded constants instead.
2235 call getenv_loc('SCPPAR',scpname)
2236 #if defined(WINIFL) || defined(WINPGI)
2237 open (iscpp,file=scpname,status='old',readonly,shared)
2238 #elif (defined CRAY) || (defined AIX)
2239 open (iscpp,file=scpname,status='old',action='read')
2241 open (iscpp,file=scpname,status='old')
2243 open (iscpp,file=scpname,status='old',readonly)
2246 call getenv_loc('PATTERN',patname)
2247 #if defined(WINIFL) || defined(WINPGI)
2248 open (icbase,file=patname,status='old',readonly,shared)
2249 #elif (defined CRAY) || (defined AIX)
2250 open (icbase,file=patname,status='old',action='read')
2252 open (icbase,file=patname,status='old')
2254 open (icbase,file=patname,status='old',readonly)
2257 C Open output file only for CG processes
2258 c print *,"Processor",myrank," fg_rank",fg_rank
2259 if (fg_rank.eq.0) then
2261 if (nodes.eq.1) then
2264 npos = dlog10(dfloat(nodes-1))+1
2266 if (npos.lt.3) npos=3
2267 write (liczba,'(i1)') npos
2268 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2270 write (liczba,form) me
2271 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2272 & liczba(:ilen(liczba))
2273 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2275 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2277 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2278 & liczba(:ilen(liczba))//'.mol2'
2279 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2280 & liczba(:ilen(liczba))//'.stat'
2282 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2283 & //liczba(:ilen(liczba))//'.stat')
2284 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2287 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2288 & liczba(:ilen(liczba))//'.const'
2293 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2294 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2295 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2296 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2297 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2299 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2301 rest2name=prefix(:ilen(prefix))//'.rst'
2303 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2306 #if defined(AIX) || defined(PGI)
2307 if (me.eq.king .or. .not. out1file)
2308 & open(iout,file=outname,status='unknown')
2310 if (fg_rank.gt.0) then
2311 write (liczba,'(i3.3)') myrank/nfgtasks
2312 write (ll,'(bz,i3.3)') fg_rank
2313 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2318 open(igeom,file=intname,status='unknown',position='append')
2319 open(ipdb,file=pdbname,status='unknown')
2320 open(imol2,file=mol2name,status='unknown')
2321 open(istat,file=statname,status='unknown',position='append')
2323 c1out open(iout,file=outname,status='unknown')
2326 if (me.eq.king .or. .not.out1file)
2327 & open(iout,file=outname,status='unknown')
2329 if (fg_rank.gt.0) then
2330 write (liczba,'(i3.3)') myrank/nfgtasks
2331 write (ll,'(bz,i3.3)') fg_rank
2332 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2337 open(igeom,file=intname,status='unknown',access='append')
2338 open(ipdb,file=pdbname,status='unknown')
2339 open(imol2,file=mol2name,status='unknown')
2340 open(istat,file=statname,status='unknown',access='append')
2342 c1out open(iout,file=outname,status='unknown')
2345 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2346 csa_seed=prefix(:lenpre)//'.CSA.seed'
2347 csa_history=prefix(:lenpre)//'.CSA.history'
2348 csa_bank=prefix(:lenpre)//'.CSA.bank'
2349 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2350 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2351 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2352 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2353 csa_int=prefix(:lenpre)//'.int'
2354 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2355 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2356 csa_in=prefix(:lenpre)//'.CSA.in'
2357 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2360 write (iout,'(80(1h-))')
2361 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2362 write (iout,'(80(1h-))')
2363 write (iout,*) "Input file : ",
2364 & pref_orig(:ilen(pref_orig))//'.inp'
2365 write (iout,*) "Output file : ",
2366 & outname(:ilen(outname))
2368 write (iout,*) "Sidechain potential file : ",
2369 & sidename(:ilen(sidename))
2371 write (iout,*) "SCp potential file : ",
2372 & scpname(:ilen(scpname))
2374 write (iout,*) "Electrostatic potential file : ",
2375 & elename(:ilen(elename))
2376 write (iout,*) "Cumulant coefficient file : ",
2377 & fouriername(:ilen(fouriername))
2378 write (iout,*) "Torsional parameter file : ",
2379 & torname(:ilen(torname))
2380 write (iout,*) "Double torsional parameter file : ",
2381 & tordname(:ilen(tordname))
2382 write (iout,*) "SCCOR parameter file : ",
2383 & sccorname(:ilen(sccorname))
2384 write (iout,*) "Bond & inertia constant file : ",
2385 & bondname(:ilen(bondname))
2386 write (iout,*) "Bending parameter file : ",
2387 & thetname(:ilen(thetname))
2388 write (iout,*) "Rotamer parameter file : ",
2389 & rotname(:ilen(rotname))
2390 write (iout,*) "Thetpdb parameter file : ",
2391 & thetname_pdb(:ilen(thetname_pdb))
2392 write (iout,*) "Threading database : ",
2393 & patname(:ilen(patname))
2395 &write (iout,*)" DIRTMP : ",
2397 write (iout,'(80(1h-))')
2401 c----------------------------------------------------------------------------
2402 subroutine card_concat(card)
2403 implicit real*8 (a-h,o-z)
2404 include 'DIMENSIONS'
2405 include 'COMMON.IOUNITS'
2407 character*80 karta,ucase
2409 read (inp,'(a)') karta
2412 do while (karta(80:80).eq.'&')
2413 card=card(:ilen(card)+1)//karta(:79)
2414 read (inp,'(a)') karta
2417 card=card(:ilen(card)+1)//karta
2420 c----------------------------------------------------------------------------------
2422 implicit real*8 (a-h,o-z)
2423 include 'DIMENSIONS'
2424 include 'COMMON.CHAIN'
2425 include 'COMMON.IOUNITS'
2427 open(irest2,file=rest2name,status='unknown')
2428 read(irest2,*) totT,EK,potE,totE,t_bath
2431 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2433 C WARING IN OUTER FIELD CRITICAL CORRECTION TO CHANGE i=0
2435 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2438 read (irest2,*) iset
2443 c---------------------------------------------------------------------------------
2444 subroutine read_fragments
2445 implicit real*8 (a-h,o-z)
2446 include 'DIMENSIONS'
2450 include 'COMMON.SETUP'
2451 include 'COMMON.CHAIN'
2452 include 'COMMON.IOUNITS'
2454 include 'COMMON.CONTROL'
2455 read(inp,*) nset,nfrag,npair,nfrag_back
2456 if(me.eq.king.or..not.out1file)
2457 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2458 & " nfrag_back",nfrag_back
2460 read(inp,*) mset(iset)
2462 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2464 if(me.eq.king.or..not.out1file)
2465 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2466 & ifrag(2,i,iset), qinfrag(i,iset)
2469 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2471 if(me.eq.king.or..not.out1file)
2472 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2473 & ipair(2,i,iset), qinpair(i,iset)
2476 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2477 & wfrag_back(3,i,iset),
2478 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2479 if(me.eq.king.or..not.out1file)
2480 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2481 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2486 C---------------------------------------------------------------------------
2487 subroutine read_afminp
2488 implicit real*8 (a-h,o-z)
2489 include 'DIMENSIONS'
2493 include 'COMMON.SETUP'
2494 include 'COMMON.CONTROL'
2495 include 'COMMON.CHAIN'
2496 include 'COMMON.IOUNITS'
2497 include 'COMMON.SBRIDGE'
2498 character*320 afmcard
2500 call card_concat(afmcard)
2501 call readi(afmcard,"BEG",afmbeg,0)
2502 call readi(afmcard,"END",afmend,0)
2503 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2504 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2505 print *,'FORCE=' ,forceAFMconst
2506 CCCC NOW PROPERTIES FOR AFM
2509 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2511 distafminit=dsqrt(distafminit)
2512 print *,'initdist',distafminit
2516 c-------------------------------------------------------------------------------
2517 subroutine read_dist_constr
2518 implicit real*8 (a-h,o-z)
2519 include 'DIMENSIONS'
2523 include 'COMMON.SETUP'
2524 include 'COMMON.CONTROL'
2525 include 'COMMON.CHAIN'
2526 include 'COMMON.IOUNITS'
2527 include 'COMMON.SBRIDGE'
2528 integer ifrag_(2,100),ipair_(2,100)
2529 double precision wfrag_(100),wpair_(100)
2530 character*500 controlcard
2532 write (iout,*) "Calling read_dist_constr"
2533 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2535 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
2536 call gen_dist_constr
2539 call card_concat(controlcard)
2540 call readi(controlcard,"NFRAG",nfrag_,0)
2541 call readi(controlcard,"NPAIR",npair_,0)
2542 call readi(controlcard,"NDIST",ndist_,0)
2543 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2544 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2545 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2546 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2547 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2548 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2549 c write (iout,*) "IFRAG"
2551 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2553 c write (iout,*) "IPAIR"
2555 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2559 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2560 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2561 & ifrag_(2,i)=nstart_sup+nsup-1
2562 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2564 if (wfrag_(i).gt.0.0d0) then
2565 do j=ifrag_(1,i),ifrag_(2,i)-1
2566 do k=j+1,ifrag_(2,i)
2567 c write (iout,*) "j",j," k",k
2569 if (constr_dist.eq.1) then
2574 forcon(nhpb)=wfrag_(i)
2575 else if (constr_dist.eq.2) then
2576 if (ddjk.le.dist_cut) then
2581 forcon(nhpb)=wfrag_(i)
2588 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2591 if (.not.out1file .or. me.eq.king)
2592 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2593 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2595 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2596 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2603 if (wpair_(i).gt.0.0d0) then
2611 do j=ifrag_(1,ii),ifrag_(2,ii)
2612 do k=ifrag_(1,jj),ifrag_(2,jj)
2616 forcon(nhpb)=wpair_(i)
2617 dhpb(nhpb)=dist(j,k)
2619 if (.not.out1file .or. me.eq.king)
2620 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2621 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2623 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2624 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2632 if (constr_dist.eq.11) then
2633 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2634 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2635 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2638 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2639 & ibecarb(i),forcon(nhpb+1)
2641 if (forcon(nhpb+1).gt.0.0d0) then
2643 if (ibecarb(i).gt.0) then
2644 ihpb(i)=ihpb(i)+nres
2645 jhpb(i)=jhpb(i)+nres
2647 if (dhpb(nhpb).eq.0.0d0)
2648 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2650 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2651 C if (forcon(nhpb+1).gt.0.0d0) then
2653 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2655 if (.not.out1file .or. me.eq.king)
2656 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2657 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2659 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2660 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2668 c-------------------------------------------------------------------------------
2670 subroutine flush(iu)
2675 subroutine flush(iu)
2680 c------------------------------------------------------------------------------
2681 subroutine copy_to_tmp(source)
2682 include "DIMENSIONS"
2683 include "COMMON.IOUNITS"
2684 character*(*) source
2685 character* 256 tmpfile
2689 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2690 inquire(file=tmpfile,exist=ex)
2692 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2693 & " to temporary directory..."
2694 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2695 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2699 c------------------------------------------------------------------------------
2700 subroutine move_from_tmp(source)
2701 include "DIMENSIONS"
2702 include "COMMON.IOUNITS"
2703 character*(*) source
2706 write (*,*) "Moving ",source(:ilen(source)),
2707 & " from temporary directory to working directory"
2708 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2709 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2712 c------------------------------------------------------------------------------
2713 subroutine random_init(seed)
2715 C Initialize random number generator
2717 implicit real*8 (a-h,o-z)
2718 include 'DIMENSIONS'
2721 logical OKRandom, prng_restart
2723 integer iseed_array(4)
2725 include 'COMMON.IOUNITS'
2726 include 'COMMON.TIME1'
2727 include 'COMMON.THREAD'
2728 include 'COMMON.SBRIDGE'
2729 include 'COMMON.CONTROL'
2730 include 'COMMON.MCM'
2731 include 'COMMON.MAP'
2732 include 'COMMON.HEADER'
2733 include 'COMMON.CSA'
2734 include 'COMMON.CHAIN'
2735 include 'COMMON.MUCA'
2737 include 'COMMON.FFIELD'
2738 include 'COMMON.SETUP'
2739 iseed=-dint(dabs(seed))
2740 if (iseed.eq.0) then
2741 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2742 & 'Random seed undefined. The program will stop.'
2743 write (*,'(/80(1h*)/20x,a/80(1h*))')
2744 & 'Random seed undefined. The program will stop.'
2746 call mpi_finalize(mpi_comm_world,ierr)
2748 stop 'Bad random seed.'
2751 if (fg_rank.eq.0) then
2755 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2756 OKRandom = prng_restart(me,iseed)
2759 tmp=65536.0d0**(4-i)
2760 iseed_array(i) = dint(seed/tmp)
2761 seed=seed-iseed_array(i)*tmp
2764 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2765 & (iseed_array(i),i=1,4)
2766 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2767 & (iseed_array(i),i=1,4)
2768 OKRandom = prng_restart(me,iseed_array)
2771 c r1 = prng_next(me)
2772 r1=ran_number(0.0D0,1.0D0)
2774 & write (iout,*) 'ran_num',r1
2775 if (r1.lt.0.0d0) OKRandom=.false.
2777 if (.not.OKRandom) then
2778 write (iout,*) 'PRNG IS NOT WORKING!!!'
2779 print *,'PRNG IS NOT WORKING!!!'
2782 call mpi_abort(mpi_comm_world,error_msg,ierr)
2785 write (iout,*) 'too many processors for parallel prng'
2786 write (*,*) 'too many processors for parallel prng'
2794 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)