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 job setup parameters
15 C Read force-field parameters except weights
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 energy-term weights and disulfide parameters
39 C Read molecule information, molecule geometry, energy-term weights, and
40 C restraints if requested
42 C Print restraint information
44 if (.not. out1file .or. me.eq.king) then
47 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
48 & "The following",nhpb-nss,
49 & " distance restraints have been imposed:",
50 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
53 write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
54 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
59 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
60 & "The following",npeak,
61 & " NMR peak restraints have been imposed:",
62 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
63 & " score"," type"," ipeak"
65 do j=ipeak(1,i),ipeak(2,i)
66 write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
67 & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
68 & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
71 write (iout,*) "The ipeak array"
73 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
79 c print *,"Processor",myrank," leaves READRTNS"
82 C-------------------------------------------------------------------------------
83 subroutine read_control
87 implicit real*8 (a-h,o-z)
91 logical OKRandom, prng_restart
94 include 'COMMON.IOUNITS'
95 include 'COMMON.TIME1'
96 include 'COMMON.THREAD'
97 include 'COMMON.SBRIDGE'
98 include 'COMMON.CONTROL'
101 include 'COMMON.HEADER'
103 include 'COMMON.CHAIN'
104 include 'COMMON.MUCA'
106 include 'COMMON.FFIELD'
107 include 'COMMON.INTERACT'
108 include 'COMMON.SETUP'
109 include 'COMMON.SPLITELE'
110 include 'COMMON.SHIELD'
112 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
113 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
115 character*320 controlcard
120 read (INP,'(a)') titel
121 call card_concat(controlcard)
122 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
123 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
124 call reada(controlcard,'SEED',seed,0.0D0)
125 call random_init(seed)
126 C Set up the time limit (caution! The time must be input in minutes!)
127 read_cart=index(controlcard,'READ_CART').gt.0
128 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
129 C this variable with_theta_constr is the variable which allow to read and execute the
130 C constrains on theta angles WITH_THETA_CONSTR is the keyword
131 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
132 write (iout,*) "with_theta_constr ",with_theta_constr
133 call readi(controlcard,'NSAXS',nsaxs,0)
134 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
135 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
136 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
137 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
138 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
139 call readi(controlcard,'SYM',symetr,1)
140 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
141 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
142 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
143 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
144 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
145 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
146 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
147 call reada(controlcard,'DRMS',drms,0.1D0)
148 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
149 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
150 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
151 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
152 write (iout,'(a,f10.1)')'DRMS = ',drms
153 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
154 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
156 call readi(controlcard,'NZ_START',nz_start,0)
157 call readi(controlcard,'NZ_END',nz_end,0)
158 c call readi(controlcard,'IZ_SC',iz_sc,0)
160 safety = 60.0d0*safety
163 call reada(controlcard,"T_BATH",t_bath,300.0d0)
164 minim=(index(controlcard,'MINIMIZE').gt.0)
165 dccart=(index(controlcard,'CART').gt.0)
166 overlapsc=(index(controlcard,'OVERLAP').gt.0)
167 overlapsc=.not.overlapsc
168 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
169 searchsc=.not.searchsc
170 sideadd=(index(controlcard,'SIDEADD').gt.0)
171 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
172 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
173 outpdb=(index(controlcard,'PDBOUT').gt.0)
174 outmol2=(index(controlcard,'MOL2OUT').gt.0)
175 pdbref=(index(controlcard,'PDBREF').gt.0)
176 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
177 indpdb=index(controlcard,'PDBSTART')
178 extconf=(index(controlcard,'EXTCONF').gt.0)
179 AFMlog=(index(controlcard,'AFM'))
180 selfguide=(index(controlcard,'SELFGUIDE'))
181 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
182 call readi(controlcard,'TUBEMOD',tubelog,0)
183 c write (iout,*) TUBElog,"TUBEMODE"
184 call readi(controlcard,'IPRINT',iprint,0)
185 C SHIELD keyword sets if the shielding effect of side-chains is used
186 C 0 denots no shielding is used all peptide are equally despite the
187 C solvent accesible area
188 C 1 the newly introduced function
189 C 2 reseved for further possible developement
190 call readi(controlcard,'SHIELD',shield_mode,0)
191 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
192 write(iout,*) "shield_mode",shield_mode
194 call readi(controlcard,'TORMODE',tor_mode,0)
195 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
196 write(iout,*) "torsional and valence angle mode",tor_mode
197 call readi(controlcard,'MAXGEN',maxgen,10000)
198 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
199 call readi(controlcard,"KDIAG",kdiag,0)
200 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
201 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
202 & write (iout,*) "RESCALE_MODE",rescale_mode
203 split_ene=index(controlcard,'SPLIT_ENE').gt.0
204 if (index(controlcard,'REGULAR').gt.0.0D0) then
205 call reada(controlcard,'WEIDIS',weidis,0.1D0)
210 if (index(controlcard,"CASC").gt.0) then
212 else if (index(controlcard,"CAONLY").gt.0) then
214 else if (index(controlcard,"SCONLY").gt.0) then
217 write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
218 & " specifying REFSTR, PDBREF or REGULAR."
222 if (index(controlcard,'CHECKGRAD').gt.0) then
224 if (index(controlcard,' CART').gt.0) then
226 elseif (index(controlcard,'CARINT').gt.0) then
231 call reada(controlcard,'DELTA',aincr,1.0d-7)
232 c write (iout,*) "icheckgrad",icheckgrad
233 elseif (index(controlcard,'THREAD').gt.0) then
235 call readi(controlcard,'THREAD',nthread,0)
236 if (nthread.gt.0) then
237 call reada(controlcard,'WEIDIS',weidis,0.1D0)
240 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
241 stop 'Error termination in Read_Control.'
243 else if (index(controlcard,'MCMA').gt.0) then
245 else if (index(controlcard,'MCEE').gt.0) then
247 else if (index(controlcard,'MULTCONF').gt.0) then
249 else if (index(controlcard,'MAP').gt.0) then
251 call readi(controlcard,'MAP',nmap,0)
252 else if (index(controlcard,'CSA').gt.0) then
254 crc else if (index(controlcard,'ZSCORE').gt.0) then
256 crc ZSCORE is rm from UNRES, modecalc=9 is available
259 cfcm else if (index(controlcard,'MCMF').gt.0) then
261 else if (index(controlcard,'SOFTREG').gt.0) then
263 else if (index(controlcard,'CHECK_BOND').gt.0) then
265 else if (index(controlcard,'TEST').gt.0) then
267 else if (index(controlcard,'MD').gt.0) then
269 else if (index(controlcard,'RE ').gt.0) then
273 lmuca=index(controlcard,'MUCA').gt.0
274 call readi(controlcard,'MUCADYN',mucadyn,0)
275 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
276 if (lmuca .and. (me.eq.king .or. .not.out1file ))
278 write (iout,*) 'MUCADYN=',mucadyn
279 write (iout,*) 'MUCASMOOTH=',muca_smooth
282 iscode=index(controlcard,'ONE_LETTER')
283 indphi=index(controlcard,'PHI')
284 indback=index(controlcard,'BACK')
285 iranconf=index(controlcard,'RAND_CONF')
286 i2ndstr=index(controlcard,'USE_SEC_PRED')
287 gradout=index(controlcard,'GRADOUT').gt.0
288 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
289 C DISTCHAINMAX become obsolete for periodic boundry condition
290 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
291 C Reading the dimensions of box in x,y,z coordinates
292 call reada(controlcard,'BOXX',boxxsize,100.0d0)
293 call reada(controlcard,'BOXY',boxysize,100.0d0)
294 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
295 c Cutoff range for interactions
296 call reada(controlcard,"R_CUT",r_cut,15.0d0)
297 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
298 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
299 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
300 if (lipthick.gt.0.0d0) then
301 bordliptop=(boxzsize+lipthick)/2.0
302 bordlipbot=bordliptop-lipthick
304 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
305 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
306 buflipbot=bordlipbot+lipbufthick
307 bufliptop=bordliptop-lipbufthick
308 if ((lipbufthick*2.0d0).gt.lipthick)
309 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
311 write(iout,*) "bordliptop=",bordliptop
312 write(iout,*) "bordlipbot=",bordlipbot
313 write(iout,*) "bufliptop=",bufliptop
314 write(iout,*) "buflipbot=",buflipbot
315 write (iout,*) "SHIELD MODE",shield_mode
316 if (TUBElog.gt.0) then
317 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
318 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
319 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
320 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
321 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
322 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
323 buftubebot=bordtubebot+tubebufthick
324 buftubetop=bordtubetop-tubebufthick
326 c if (shield_mode.gt.0) then
328 C VSolvSphere the volume of solving sphere
330 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
331 C there will be no distinction between proline peptide group and normal peptide
332 C group in case of shielding parameters
333 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
334 c VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
335 c VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
336 c write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
338 C long axis of side chain
340 c long_r_sidechain(i)=vbldsc0(1,i)
341 c short_r_sidechain(i)=sigma0(i)
345 if (me.eq.king .or. .not.out1file )
346 & write (iout,*) "DISTCHAINMAX",distchainmax
348 if(me.eq.king.or..not.out1file)
349 & write (iout,'(2a)') diagmeth(kdiag),
350 & ' routine used to diagonalize matrices.'
353 c--------------------------------------------------------------------------
354 subroutine read_REMDpar
358 implicit real*8 (a-h,o-z)
360 include 'COMMON.IOUNITS'
361 include 'COMMON.TIME1'
364 include 'COMMON.LANGEVIN'
366 include 'COMMON.LANGEVIN.lang0'
368 include 'COMMON.INTERACT'
369 include 'COMMON.NAMES'
371 include 'COMMON.REMD'
372 include 'COMMON.CONTROL'
373 include 'COMMON.SETUP'
375 character*320 controlcard
376 character*3200 controlcard1
377 integer iremd_m_total
379 if(me.eq.king.or..not.out1file)
380 & write (iout,*) "REMD setup"
382 call card_concat(controlcard)
383 call readi(controlcard,"NREP",nrep,3)
384 call readi(controlcard,"NSTEX",nstex,1000)
385 call reada(controlcard,"RETMIN",retmin,10.0d0)
386 call reada(controlcard,"RETMAX",retmax,1000.0d0)
387 mremdsync=(index(controlcard,'SYNC').gt.0)
388 call readi(controlcard,"NSYN",i_sync_step,100)
389 restart1file=(index(controlcard,'REST1FILE').gt.0)
390 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
391 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
392 if(max_cache_traj_use.gt.max_cache_traj)
393 & max_cache_traj_use=max_cache_traj
394 if(me.eq.king.or..not.out1file) then
395 cd if (traj1file) then
396 crc caching is in testing - NTWX is not ignored
397 cd write (iout,*) "NTWX value is ignored"
398 cd write (iout,*) " trajectory is stored to one file by master"
399 cd write (iout,*) " before exchange at NSTEX intervals"
401 write (iout,*) "NREP= ",nrep
402 write (iout,*) "NSTEX= ",nstex
403 write (iout,*) "SYNC= ",mremdsync
404 write (iout,*) "NSYN= ",i_sync_step
405 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
408 if (index(controlcard,'TLIST').gt.0) then
410 call card_concat(controlcard1)
411 read(controlcard1,*) (remd_t(i),i=1,nrep)
412 if(me.eq.king.or..not.out1file)
413 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
416 if (index(controlcard,'MLIST').gt.0) then
418 call card_concat(controlcard1)
419 read(controlcard1,*) (remd_m(i),i=1,nrep)
420 if(me.eq.king.or..not.out1file) then
421 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
424 iremd_m_total=iremd_m_total+remd_m(i)
426 write (iout,*) 'Total number of replicas ',iremd_m_total
431 & "Adaptive (PMF-biased) umbrella sampling will be run"
434 if(me.eq.king.or..not.out1file)
435 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
438 c--------------------------------------------------------------------------
439 subroutine read_MDpar
443 implicit real*8 (a-h,o-z)
445 include 'COMMON.IOUNITS'
446 include 'COMMON.TIME1'
449 include 'COMMON.LANGEVIN'
451 include 'COMMON.LANGEVIN.lang0'
453 include 'COMMON.INTERACT'
454 include 'COMMON.NAMES'
456 include 'COMMON.SETUP'
457 include 'COMMON.CONTROL'
458 include 'COMMON.SPLITELE'
459 include 'COMMON.FFIELD'
461 character*320 controlcard
463 call card_concat(controlcard)
464 call readi(controlcard,"NSTEP",n_timestep,1000000)
465 call readi(controlcard,"NTWE",ntwe,100)
466 call readi(controlcard,"NTWX",ntwx,1000)
467 call reada(controlcard,"DT",d_time,1.0d-1)
468 call reada(controlcard,"DVMAX",dvmax,2.0d1)
469 call reada(controlcard,"DAMAX",damax,1.0d1)
470 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
471 call readi(controlcard,"LANG",lang,0)
472 RESPA = index(controlcard,"RESPA") .gt. 0
473 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
474 ntime_split0=ntime_split
475 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
476 ntime_split0=ntime_split
477 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
478 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
479 rest = index(controlcard,"REST").gt.0
480 tbf = index(controlcard,"TBF").gt.0
481 usampl = index(controlcard,"USAMPL").gt.0
482 scale_umb = index(controlcard,"SCALE_UMB").gt.0
483 adaptive = index(controlcard,"ADAPTIVE").gt.0
484 mdpdb = index(controlcard,"MDPDB").gt.0
485 call reada(controlcard,"T_BATH",t_bath,300.0d0)
486 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
487 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
488 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
489 if (count_reset_moment.eq.0) count_reset_moment=1000000000
490 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
491 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
492 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
493 if (count_reset_vel.eq.0) count_reset_vel=1000000000
494 large = index(controlcard,"LARGE").gt.0
495 print_compon = index(controlcard,"PRINT_COMPON").gt.0
496 rattle = index(controlcard,"RATTLE").gt.0
497 c if performing umbrella sampling, fragments constrained are read from the fragment file
500 write (iout,*) "Umbrella sampling will be run"
501 if (scale_umb.and.adaptive) then
502 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
503 write (iout,*) "Select one of those and re-run the job."
506 if (scale_umb) write (iout,*)
507 &"Umbrella-restraint force constants will be scaled by temperature"
511 if(me.eq.king.or..not.out1file) then
513 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
515 write (iout,'(a)') "The units are:"
516 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
517 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
518 & " acceleration: angstrom/(48.9 fs)**2"
519 write (iout,'(a)') "energy: kcal/mol, temperature: K"
521 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
522 write (iout,'(a60,f10.5,a)')
523 & "Initial time step of numerical integration:",d_time,
525 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
527 write (iout,'(2a,i4,a)')
528 & "A-MTS algorithm used; initial time step for fast-varying",
529 & " short-range forces split into",ntime_split," steps."
530 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
531 & r_cut," lambda",rlamb
533 write (iout,'(2a,f10.5)')
534 & "Maximum acceleration threshold to reduce the time step",
535 & "/increase split number:",damax
536 write (iout,'(2a,f10.5)')
537 & "Maximum predicted energy drift to reduce the timestep",
538 & "/increase split number:",edriftmax
539 write (iout,'(a60,f10.5)')
540 & "Maximum velocity threshold to reduce velocities:",dvmax
541 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
542 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
543 if (rattle) write (iout,'(a60)')
544 & "Rattle algorithm used to constrain the virtual bonds"
548 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
549 call reada(controlcard,"RWAT",rwat,1.4d0)
550 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
551 surfarea=index(controlcard,"SURFAREA").gt.0
552 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
553 if(me.eq.king.or..not.out1file)then
554 write (iout,'(/a,$)') "Langevin dynamics calculation"
557 & " with direct integration of Langevin equations"
558 else if (lang.eq.2) then
559 write (iout,'(a/)') " with TINKER stochasic MD integrator"
560 else if (lang.eq.3) then
561 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
562 else if (lang.eq.4) then
563 write (iout,'(a/)') " in overdamped mode"
565 write (iout,'(//a,i5)')
566 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
569 write (iout,'(a60,f10.5)') "Temperature:",t_bath
570 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
571 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
572 write (iout,'(a60,f10.5)')
573 & "Scaling factor of the friction forces:",scal_fric
574 if (surfarea) write (iout,'(2a,i10,a)')
575 & "Friction coefficients will be scaled by solvent-accessible",
576 & " surface area every",reset_fricmat," steps."
578 c Calculate friction coefficients and bounds of stochastic forces
579 eta=6*pi*cPoise*etawat
580 if(me.eq.king.or..not.out1file)
581 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
583 gamp=scal_fric*(pstok+rwat)*eta
584 stdfp=dsqrt(2*Rb*t_bath/d_time)
586 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
587 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
589 if(me.eq.king.or..not.out1file)then
590 write (iout,'(/2a/)')
591 & "Radii of site types and friction coefficients and std's of",
592 & " stochastic forces of fully exposed sites"
593 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
595 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
596 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
600 if(me.eq.king.or..not.out1file)then
601 write (iout,'(a)') "Berendsen bath calculation"
602 write (iout,'(a60,f10.5)') "Temperature:",t_bath
603 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
605 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
606 & count_reset_moment," steps"
608 & write (iout,'(a,i10,a)')
609 & "Velocities will be reset at random every",count_reset_vel,
613 if(me.eq.king.or..not.out1file)
614 & write (iout,'(a31)') "Microcanonical mode calculation"
616 if(me.eq.king.or..not.out1file)then
617 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
619 write(iout,*) "MD running with constraints."
620 write(iout,*) "Equilibration time ", eq_time, " mtus."
621 write(iout,*) "Constraining ", nfrag," fragments."
622 write(iout,*) "Length of each fragment, weight and q0:"
624 write (iout,*) "Set of restraints #",iset
626 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
627 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
629 write(iout,*) "constraints between ", npair, "fragments."
630 write(iout,*) "constraint pairs, weights and q0:"
632 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
633 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
635 write(iout,*) "angle constraints within ", nfrag_back,
636 & "backbone fragments."
638 write(iout,*) "fragment, weights, q0:"
640 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
641 & ifrag_back(2,i,iset),
642 & wfrag_back(1,i,iset),qin_back(1,i,iset),
643 & wfrag_back(2,i,iset),qin_back(2,i,iset),
644 & wfrag_back(3,i,iset),qin_back(3,i,iset)
647 write(iout,*) "fragment, weights:"
649 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
650 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
651 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
655 iset=mod(kolor,nset)+1
658 if(me.eq.king.or..not.out1file)
659 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
662 c------------------------------------------------------------------------------
665 C Read molecular data.
667 implicit real*8 (a-h,o-z)
673 include 'COMMON.IOUNITS'
676 include 'COMMON.INTERACT'
677 include 'COMMON.LOCAL'
678 include 'COMMON.NAMES'
679 include 'COMMON.CHAIN'
680 include 'COMMON.FFIELD'
681 include 'COMMON.SBRIDGE'
682 include 'COMMON.HEADER'
683 include 'COMMON.CONTROL'
684 include 'COMMON.DBASE'
685 include 'COMMON.THREAD'
686 include 'COMMON.CONTACTS'
687 include 'COMMON.TORCNSTR'
688 include 'COMMON.TIME1'
689 include 'COMMON.BOUNDS'
691 include 'COMMON.SETUP'
692 include 'COMMON.SHIELD'
693 character*4 sequence(maxres)
695 double precision x(maxvar)
696 character*256 pdbfile
697 character*400 weightcard
698 character*80 weightcard_t,ucase
699 dimension itype_pdb(maxres)
700 common /pizda/ itype_pdb
701 logical seq_comp,fail
702 double precision energia(0:n_ene)
703 double precision secprob(3,maxdih_constr)
707 C Read PDB structure if applicable
709 if (indpdb.gt.0 .or. pdbref) then
710 read(inp,'(a)') pdbfile
711 if(me.eq.king.or..not.out1file)
712 & write (iout,'(2a)') 'PDB data will be read from file ',
713 & pdbfile(:ilen(pdbfile))
714 open(ipdbin,file=pdbfile,status='old',err=33)
716 33 write (iout,'(a)') 'Error opening PDB file.'
720 if(me.eq.king.or..not.out1file)
721 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
722 & ' nstart_sup=',nstart_sup
724 itype_pdb(i)=itype(i)
728 nct=nstart_sup+nsup-1
729 call contact(.false.,ncont_ref,icont_ref,co)
732 if(me.eq.king.or..not.out1file)
733 & write(iout,*)'Adding sidechains'
737 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
740 do while (fail.and.nsi.le.maxsi)
741 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
744 if(fail) write(iout,*)'Adding sidechain failed for res ',
745 & i,' after ',nsi,' trials'
750 if (indpdb.eq.0) then
751 C Read sequence if not taken from the pdb file.
753 c print *,'nres=',nres
754 if (iscode.gt.0) then
755 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
757 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
759 C Convert sequence to numeric code
761 itype(i)=rescode(i,sequence(i),iscode)
763 C Assign initial virtual bond lengths
769 vbld(i+nres)=dsc(iabs(itype(i)))
770 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
771 c write (iout,*) "i",i," itype",itype(i),
772 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
776 c print '(20i4)',(itype(i),i=1,nres)
779 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
781 if (itype(i).eq.ntyp1) then
785 else if (iabs(itype(i+1)).ne.20) then
787 else if (iabs(itype(i)).ne.20) then
794 if(me.eq.king.or..not.out1file)then
795 write (iout,*) "ITEL"
797 write (iout,*) i,itype(i),itel(i)
799 print *,'Call Read_Bridge.'
803 cd print *,'NNT=',NNT,' NCT=',NCT
804 call seq2chains(nres,itype,nchain,chain_length,chain_border,
806 write(iout,*) "nres",nres," nchain",nchain
808 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
811 call chain_symmetry(nchain,nres,itype,chain_border,
812 & chain_length,npermchain,tabpermchain)
814 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
817 write(iout,*) "residue permutations"
819 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
821 if (itype(1).eq.ntyp1) nnt=2
822 if (itype(nres).eq.ntyp1) nct=nct-1
824 if(me.eq.king.or..not.out1file)
825 & write (iout,'(a,i3)') 'nsup=',nsup
827 if (nsup.le.(nct-nnt+1)) then
828 do i=0,nct-nnt+1-nsup
829 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
835 & 'Error - sequences to be superposed do not match.'
838 do i=0,nsup-(nct-nnt+1)
839 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
841 nstart_sup=nstart_sup+i
847 & 'Error - sequences to be superposed do not match.'
850 if (nsup.eq.0) nsup=nct-nnt
851 if (nstart_sup.eq.0) nstart_sup=nnt
852 if (nstart_seq.eq.0) nstart_seq=nnt
853 if(me.eq.king.or..not.out1file)
854 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
855 & ' nstart_seq=',nstart_seq
858 C 8/13/98 Set limits to generating the dihedral angles
863 read (inp,*) ndih_constr
864 write (iout,*) "ndih_constr",ndih_constr
865 if (ndih_constr.gt.0) then
868 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
871 if(me.eq.king.or..not.out1file)then
874 & 'There are',ndih_constr,' restraints on gamma angles.'
876 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
883 phi0(i)=deg2rad*phi0(i)
884 drange(i)=deg2rad*drange(i)
886 C if(me.eq.king.or..not.out1file)
887 C & write (iout,*) 'FTORS',ftors
890 phibound(1,ii) = phi0(i)-drange(i)
891 phibound(2,ii) = phi0(i)+drange(i)
894 if (me.eq.king .or. .not.out1file) then
896 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
898 write (iout,'(a3,i5,2f10.1)')
899 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
904 else if (ndih_constr.lt.0) then
906 call card_concat(weightcard)
907 call reada(weightcard,"PHIHEL",phihel,50.0D0)
908 call reada(weightcard,"PHIBET",phibet,180.0D0)
909 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
910 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
911 call reada(weightcard,"WDIHC",wdihc,0.591D0)
912 write (iout,*) "Weight of dihedral angle restraints",wdihc
913 read(inp,'(9x,3f7.3)')
914 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
915 write (iout,*) "The secprob array"
917 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
921 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
922 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
923 ndih_constr=ndih_constr+1
924 idih_constr(ndih_constr)=i
927 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
928 sumv=sumv+vpsipred(j,ndih_constr)
931 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
933 phibound(1,ndih_constr)=phihel*deg2rad
934 phibound(2,ndih_constr)=phibet*deg2rad
935 sdihed(1,ndih_constr)=sigmahel*deg2rad
936 sdihed(2,ndih_constr)=sigmabet*deg2rad
940 if(me.eq.king.or..not.out1file)then
943 & 'There are',ndih_constr,
944 & ' bimodal restraints on gamma angles.'
946 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
947 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
948 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
949 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
950 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
951 & (vpsipred(j,i),j=1,3)
958 C first setting the theta boundaries to 0 to pi
959 C this mean that there is no energy penalty for any angle occuring this can be applied
960 C for generate random conformation but is not implemented in this way
965 C begin reading theta constrains this is quartic constrains allowing to
966 C have smooth second derivative
967 if (with_theta_constr) then
968 C with_theta_constr is keyword allowing for occurance of theta constrains
969 read (inp,*) ntheta_constr
970 C ntheta_constr is the number of theta constrains
971 if (ntheta_constr.gt.0) then
973 read (inp,*) (itheta_constr(i),theta_constr0(i),
974 & theta_drange(i),for_thet_constr(i),
976 C the above code reads from 1 to ntheta_constr
977 C itheta_constr(i) residue i for which is theta_constr
978 C theta_constr0 the global minimum value
979 C theta_drange is range for which there is no energy penalty
980 C for_thet_constr is the force constant for quartic energy penalty
982 if(me.eq.king.or..not.out1file)then
984 & 'There are',ntheta_constr,' constraints on phi angles.'
986 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
992 theta_constr0(i)=deg2rad*theta_constr0(i)
993 theta_drange(i)=deg2rad*theta_drange(i)
995 C if(me.eq.king.or..not.out1file)
996 C & write (iout,*) 'FTORS',ftors
997 C do i=1,ntheta_constr
998 C ii = itheta_constr(i)
999 C thetabound(1,ii) = phi0(i)-drange(i)
1000 C thetabound(2,ii) = phi0(i)+drange(i)
1002 endif ! ntheta_constr.gt.0
1003 endif! with_theta_constr
1005 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1006 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1007 c--- Zscore rms -------
1008 if (nz_start.eq.0) nz_start=nnt
1009 if (nz_end.eq.0 .and. nsup.gt.0) then
1011 else if (nz_end.eq.0) then
1014 if(me.eq.king.or..not.out1file)then
1015 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1016 write (iout,*) 'IZ_SC=',iz_sc
1018 c----------------------
1021 if (.not.pdbref) then
1022 call read_angles(inp,*38)
1024 38 write (iout,'(a)') 'Error reading reference structure.'
1026 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1027 stop 'Error reading reference structure'
1029 39 call chainbuild_extconf
1031 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1040 call contact(.true.,ncont_ref,icont_ref,co)
1044 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1046 if (constr_dist.gt.0) call read_dist_constr
1047 write (iout,*) "After read_dist_constr nhpb",nhpb
1048 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1050 call NMRpeak_partition
1051 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1053 if(me.eq.king.or..not.out1file)
1054 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1057 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1059 if(me.eq.king.or..not.out1file)
1060 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1061 & icont_ref(1,i),' ',
1062 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1065 write (iout,*) "calling read_saxs_consrtr",nsaxs
1066 if (nsaxs.gt.0) call read_saxs_constr
1069 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1070 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1071 & modecalc.ne.10) then
1072 C If input structure hasn't been supplied from the PDB file read or generate
1074 if (iranconf.eq.0 .and. .not. extconf) then
1075 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1076 & write (iout,'(a)') 'Initial geometry will be read in.'
1078 read(inp,'(8f10.5)',end=36,err=36)
1079 & ((c(l,k),l=1,3),k=1,nres),
1080 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1081 write (iout,*) "Exit READ_CART"
1082 c write (iout,'(8f10.5)')
1083 c & ((c(l,k),l=1,3),k=1,nres),
1084 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1088 dc(j,i)=c(j,i+1)-c(j,i)
1089 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1093 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1095 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1096 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1101 c dc_norm(j,i+nres)=0.0d0
1105 call int_from_cart1(.true.)
1106 write (iout,*) "Finish INT_TO_CART"
1107 c write (iout,*) "DC"
1109 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1110 c & (dc(j,i+nres),j=1,3)
1116 call read_angles(inp,*36)
1117 call chainbuild_extconf
1120 36 write (iout,'(a)') 'Error reading angle file.'
1122 call mpi_finalize( MPI_COMM_WORLD,IERR )
1124 stop 'Error reading angle file.'
1126 else if (extconf) then
1127 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1128 & write (iout,'(a)') 'Extended chain initial geometry.'
1130 theta(i)=90d0*deg2rad
1133 phi(i)=180d0*deg2rad
1136 alph(i)=110d0*deg2rad
1139 omeg(i)=-120d0*deg2rad
1140 if (itype(i).le.0) omeg(i)=-omeg(i)
1142 call chainbuild_extconf
1144 if(me.eq.king.or..not.out1file)
1145 & write (iout,'(a)') 'Random-generated initial geometry.'
1149 if (me.eq.king .or. fg_rank.eq.0 .and. (
1150 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1154 call gen_rand_conf(itmp,*30)
1156 30 write (iout,*) 'Failed to generate random conformation',
1157 & ', itrial=',itrial
1158 write (*,*) 'Processor:',me,
1159 & ' Failed to generate random conformation',
1169 write (iout,'(a,i3,a)') 'Processor:',me,
1170 & ' error in generating random conformation.'
1171 write (*,'(a,i3,a)') 'Processor:',me,
1172 & ' error in generating random conformation.'
1175 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1180 & ' error in generating random conformation.'
1185 elseif (modecalc.eq.4) then
1186 read (inp,'(a)') intinname
1187 open (intin,file=intinname,status='old',err=333)
1188 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1189 & write (iout,'(a)') 'intinname',intinname
1190 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1192 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1194 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1196 stop 'Error opening angle file.'
1200 C Generate distance constraints, if the PDB structure is to be regularized.
1201 if (nthread.gt.0) then
1202 call read_threadbase
1205 if (me.eq.king .or. .not. out1file)
1207 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1208 write (iout,'(/a,i3,a)')
1209 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1210 write (iout,'(20i4)') (iss(i),i=1,ns)
1212 write(iout,*)"Running with dynamic disulfide-bond formation"
1214 write (iout,'(/a/)') 'Pre-formed links are:'
1220 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1221 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1227 if (ns.gt.0.and.dyn_ss) then
1231 forcon(i-nss)=forcon(i)
1238 dyn_ss_mask(iss(i))=.true.
1243 c write (iout,*) "DC"
1245 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1246 c & (dc(j,i+nres),j=1,3)
1248 if (i2ndstr.gt.0) call secstrp2dihc
1249 c call geom_to_var(nvar,x)
1250 c call etotal(energia(0))
1251 c call enerprint(energia(0))
1252 c call briefout(0,etot)
1254 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1255 cd write (iout,'(a)') 'Variable list:'
1256 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1258 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1259 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1260 & 'Processor',myrank,': end reading molecular data.'
1265 c--------------------------------------------------------------------------
1266 logical function seq_comp(itypea,itypeb,length)
1268 integer length,itypea(length),itypeb(length)
1271 if (itypea(i).ne.itypeb(i)) then
1279 c-----------------------------------------------------------------------------
1280 subroutine read_bridge
1281 C Read information about disulfide bridges.
1282 implicit real*8 (a-h,o-z)
1283 include 'DIMENSIONS'
1287 include 'COMMON.IOUNITS'
1288 include 'COMMON.GEO'
1289 include 'COMMON.VAR'
1290 include 'COMMON.INTERACT'
1291 include 'COMMON.LOCAL'
1292 include 'COMMON.NAMES'
1293 include 'COMMON.CHAIN'
1294 include 'COMMON.FFIELD'
1295 include 'COMMON.SBRIDGE'
1296 include 'COMMON.HEADER'
1297 include 'COMMON.CONTROL'
1298 include 'COMMON.DBASE'
1299 include 'COMMON.THREAD'
1300 include 'COMMON.TIME1'
1301 include 'COMMON.SETUP'
1302 C Read bridging residues.
1303 read (inp,*) ns,(iss(i),i=1,ns)
1305 if(me.eq.king.or..not.out1file)
1306 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1307 C Check whether the specified bridging residues are cystines.
1309 if (itype(iss(i)).ne.1) then
1310 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1311 & 'Do you REALLY think that the residue ',
1312 & restyp(itype(iss(i))),i,
1313 & ' can form a disulfide bridge?!!!'
1314 write (*,'(2a,i3,a)')
1315 & 'Do you REALLY think that the residue ',
1316 & restyp(itype(iss(i))),i,
1317 & ' can form a disulfide bridge?!!!'
1319 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1324 C Read preformed bridges.
1326 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1328 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1331 C Check if the residues involved in bridges are in the specified list of
1332 C bridging residues.
1335 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1336 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1337 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1338 & ' contains residues present in other pairs.'
1339 write (*,'(a,i3,a)') 'Disulfide pair',i,
1340 & ' contains residues present in other pairs.'
1342 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1348 if (ihpb(i).eq.iss(j)) goto 10
1350 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1353 if (jhpb(i).eq.iss(j)) goto 20
1355 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1361 ihpb(i)=ihpb(i)+nres
1362 jhpb(i)=jhpb(i)+nres
1368 c----------------------------------------------------------------------------
1369 subroutine read_x(kanal,*)
1370 implicit real*8 (a-h,o-z)
1371 include 'DIMENSIONS'
1372 include 'COMMON.GEO'
1373 include 'COMMON.VAR'
1374 include 'COMMON.CHAIN'
1375 include 'COMMON.IOUNITS'
1376 include 'COMMON.CONTROL'
1377 include 'COMMON.LOCAL'
1378 include 'COMMON.INTERACT'
1379 c Read coordinates from input
1381 read(kanal,'(8f10.5)',end=10,err=10)
1382 & ((c(l,k),l=1,3),k=1,nres),
1383 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1386 c(j,2*nres)=c(j,nres)
1388 call int_from_cart1(.false.)
1391 dc(j,i)=c(j,i+1)-c(j,i)
1392 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1396 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1398 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1399 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1407 c----------------------------------------------------------------------------
1408 subroutine read_threadbase
1409 implicit real*8 (a-h,o-z)
1410 include 'DIMENSIONS'
1411 include 'COMMON.IOUNITS'
1412 include 'COMMON.GEO'
1413 include 'COMMON.VAR'
1414 include 'COMMON.INTERACT'
1415 include 'COMMON.LOCAL'
1416 include 'COMMON.NAMES'
1417 include 'COMMON.CHAIN'
1418 include 'COMMON.FFIELD'
1419 include 'COMMON.SBRIDGE'
1420 include 'COMMON.HEADER'
1421 include 'COMMON.CONTROL'
1422 include 'COMMON.DBASE'
1423 include 'COMMON.THREAD'
1424 include 'COMMON.TIME1'
1425 C Read pattern database for threading.
1426 read (icbase,*) nseq
1428 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1429 & nres_base(2,i),nres_base(3,i)
1430 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1432 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1433 c & nres_base(2,i),nres_base(3,i)
1434 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1438 if (weidis.eq.0.0D0) weidis=0.1D0
1447 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1448 write (iout,'(a,i5)') 'nexcl: ',nexcl
1449 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1452 c------------------------------------------------------------------------------
1453 subroutine setup_var
1454 implicit real*8 (a-h,o-z)
1455 include 'DIMENSIONS'
1456 include 'COMMON.IOUNITS'
1457 include 'COMMON.GEO'
1458 include 'COMMON.VAR'
1459 include 'COMMON.INTERACT'
1460 include 'COMMON.LOCAL'
1461 include 'COMMON.NAMES'
1462 include 'COMMON.CHAIN'
1463 include 'COMMON.FFIELD'
1464 include 'COMMON.SBRIDGE'
1465 include 'COMMON.HEADER'
1466 include 'COMMON.CONTROL'
1467 include 'COMMON.DBASE'
1468 include 'COMMON.THREAD'
1469 include 'COMMON.TIME1'
1470 C Set up variable list.
1475 write (iout,*) "SETUP_VAR ialph"
1477 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1479 ialph(i,1)=nvar+nside
1483 if (indphi.gt.0) then
1485 else if (indback.gt.0) then
1490 write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1493 c----------------------------------------------------------------------------
1494 subroutine gen_dist_constr
1495 C Generate CA distance constraints.
1496 implicit real*8 (a-h,o-z)
1497 include 'DIMENSIONS'
1498 include 'COMMON.IOUNITS'
1499 include 'COMMON.GEO'
1500 include 'COMMON.VAR'
1501 include 'COMMON.INTERACT'
1502 include 'COMMON.LOCAL'
1503 include 'COMMON.NAMES'
1504 include 'COMMON.CHAIN'
1505 include 'COMMON.FFIELD'
1506 include 'COMMON.SBRIDGE'
1507 include 'COMMON.HEADER'
1508 include 'COMMON.CONTROL'
1509 include 'COMMON.DBASE'
1510 include 'COMMON.THREAD'
1511 include 'COMMON.TIME1'
1512 dimension itype_pdb(maxres)
1513 common /pizda/ itype_pdb
1515 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1516 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1517 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1519 do i=nstart_sup,nstart_sup+nsup-1
1520 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1521 cd & ' seq_pdb', restyp(itype_pdb(i))
1522 do j=i+2,nstart_sup+nsup-1
1524 ihpb(nhpb)=i+nstart_seq-nstart_sup
1525 jhpb(nhpb)=j+nstart_seq-nstart_sup
1527 dhpb(nhpb)=dist(i,j)
1530 cd write (iout,'(a)') 'Distance constraints:'
1535 cd if (ii.gt.nres) then
1540 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1541 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1542 cd & dhpb(i),forcon(i)
1546 c----------------------------------------------------------------------------
1548 implicit real*8 (a-h,o-z)
1549 include 'DIMENSIONS'
1550 include 'COMMON.MAP'
1551 include 'COMMON.IOUNITS'
1552 character*3 angid(4) /'THE','PHI','ALP','OME'/
1553 character*80 mapcard,ucase
1555 read (inp,'(a)') mapcard
1556 mapcard=ucase(mapcard)
1557 if (index(mapcard,'PHI').gt.0) then
1559 else if (index(mapcard,'THE').gt.0) then
1561 else if (index(mapcard,'ALP').gt.0) then
1563 else if (index(mapcard,'OME').gt.0) then
1566 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1567 stop 'Error - illegal variable spec in MAP card.'
1569 call readi (mapcard,'RES1',res1(imap),0)
1570 call readi (mapcard,'RES2',res2(imap),0)
1571 if (res1(imap).eq.0) then
1572 res1(imap)=res2(imap)
1573 else if (res2(imap).eq.0) then
1574 res2(imap)=res1(imap)
1576 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1578 & 'Error - illegal definition of variable group in MAP.'
1579 stop 'Error - illegal definition of variable group in MAP.'
1581 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1582 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1583 call readi(mapcard,'NSTEP',nstep(imap),0)
1584 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1586 & 'Illegal boundary and/or step size specification in MAP.'
1587 stop 'Illegal boundary and/or step size specification in MAP.'
1592 c----------------------------------------------------------------------------
1594 implicit real*8 (a-h,o-z)
1595 include 'DIMENSIONS'
1596 include 'COMMON.IOUNITS'
1597 include 'COMMON.GEO'
1598 include 'COMMON.CSA'
1599 include 'COMMON.BANK'
1600 include 'COMMON.CONTROL'
1602 character*620 mcmcard
1603 call card_concat(mcmcard)
1605 call readi(mcmcard,'NCONF',nconf,50)
1606 call readi(mcmcard,'NADD',nadd,0)
1607 call readi(mcmcard,'JSTART',jstart,1)
1608 call readi(mcmcard,'JEND',jend,1)
1609 call readi(mcmcard,'NSTMAX',nstmax,500000)
1610 call readi(mcmcard,'N0',n0,1)
1611 call readi(mcmcard,'N1',n1,6)
1612 call readi(mcmcard,'N2',n2,4)
1613 call readi(mcmcard,'N3',n3,0)
1614 call readi(mcmcard,'N4',n4,0)
1615 call readi(mcmcard,'N5',n5,0)
1616 call readi(mcmcard,'N6',n6,10)
1617 call readi(mcmcard,'N7',n7,0)
1618 call readi(mcmcard,'N8',n8,0)
1619 call readi(mcmcard,'N9',n9,0)
1620 call readi(mcmcard,'N14',n14,0)
1621 call readi(mcmcard,'N15',n15,0)
1622 call readi(mcmcard,'N16',n16,0)
1623 call readi(mcmcard,'N17',n17,0)
1624 call readi(mcmcard,'N18',n18,0)
1626 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1628 call readi(mcmcard,'NDIFF',ndiff,2)
1629 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1630 call readi(mcmcard,'IS1',is1,1)
1631 call readi(mcmcard,'IS2',is2,8)
1632 call readi(mcmcard,'NRAN0',nran0,4)
1633 call readi(mcmcard,'NRAN1',nran1,2)
1634 call readi(mcmcard,'IRR',irr,1)
1635 call readi(mcmcard,'NSEED',nseed,20)
1636 call readi(mcmcard,'NTOTAL',ntotal,10000)
1637 call reada(mcmcard,'CUT1',cut1,2.0d0)
1638 call reada(mcmcard,'CUT2',cut2,5.0d0)
1639 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1640 call readi(mcmcard,'ICMAX',icmax,3)
1641 call readi(mcmcard,'IRESTART',irestart,0)
1642 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1645 call reada(mcmcard,'DELE',dele,20.0d0)
1646 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1647 call readi(mcmcard,'IREF',iref,0)
1648 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1649 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1650 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1651 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1652 write (iout,*) "NCONF_IN",nconf_in
1655 c----------------------------------------------------------------------------
1656 cfmc subroutine mcmfread
1657 cfmc implicit real*8 (a-h,o-z)
1658 cfmc include 'DIMENSIONS'
1659 cfmc include 'COMMON.MCMF'
1660 cfmc include 'COMMON.IOUNITS'
1661 cfmc include 'COMMON.GEO'
1662 cfmc character*80 ucase
1663 cfmc character*620 mcmcard
1664 cfmc call card_concat(mcmcard)
1666 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1667 cfmc write(iout,*)'MAXRANT=',maxrant
1668 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1669 cfmc write(iout,*)'MAXFAM=',maxfam
1670 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1671 cfmc write(iout,*)'NNET1=',nnet1
1672 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1673 cfmc write(iout,*)'NNET2=',nnet2
1674 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1675 cfmc write(iout,*)'NNET3=',nnet3
1676 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1677 cfmc write(iout,*)'ILASTT=',ilastt
1678 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1679 cfmc write(iout,*)'MAXSTR=',maxstr
1680 cfmc maxstr_f=maxstr/maxfam
1681 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1682 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1683 cfmc write(iout,*)'NMCMF=',nmcmf
1684 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1685 cfmc write(iout,*)'IFOCUS=',ifocus
1686 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1687 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1688 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1689 cfmc write(iout,*)'INTPRT=',intprt
1690 cfmc call readi(mcmcard,'IPRT',iprt,100)
1691 cfmc write(iout,*)'IPRT=',iprt
1692 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1693 cfmc write(iout,*)'IMAXTR=',imaxtr
1694 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1695 cfmc write(iout,*)'MAXEVEN=',maxeven
1696 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1697 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1698 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1699 cfmc write(iout,*)'INIMIN=',inimin
1700 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1701 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1702 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1703 cfmc write(iout,*)'NTHREAD=',nthread
1704 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1705 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1706 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1707 cfmc write(iout,*)'MAXPERT=',maxpert
1708 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1709 cfmc write(iout,*)'IRMSD=',irmsd
1710 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1711 cfmc write(iout,*)'DENEMIN=',denemin
1712 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1713 cfmc write(iout,*)'RCUT1S=',rcut1s
1714 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1715 cfmc write(iout,*)'RCUT1E=',rcut1e
1716 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1717 cfmc write(iout,*)'RCUT2S=',rcut2s
1718 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1719 cfmc write(iout,*)'RCUT2E=',rcut2e
1720 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1721 cfmc write(iout,*)'DPERT1=',d_pert1
1722 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1723 cfmc write(iout,*)'DPERT1A=',d_pert1a
1724 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1725 cfmc write(iout,*)'DPERT2=',d_pert2
1726 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1727 cfmc write(iout,*)'DPERT2A=',d_pert2a
1728 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1729 cfmc write(iout,*)'DPERT2B=',d_pert2b
1730 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1731 cfmc write(iout,*)'DPERT2C=',d_pert2c
1732 cfmc d_pert1=deg2rad*d_pert1
1733 cfmc d_pert1a=deg2rad*d_pert1a
1734 cfmc d_pert2=deg2rad*d_pert2
1735 cfmc d_pert2a=deg2rad*d_pert2a
1736 cfmc d_pert2b=deg2rad*d_pert2b
1737 cfmc d_pert2c=deg2rad*d_pert2c
1738 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1739 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1740 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1741 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1742 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1743 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1744 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1745 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1746 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1747 cfmc write(iout,*)'RCUTINI=',rcutini
1748 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1749 cfmc write(iout,*)'GRAT=',grat
1750 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1751 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1755 c----------------------------------------------------------------------------
1757 implicit real*8 (a-h,o-z)
1758 include 'DIMENSIONS'
1759 include 'COMMON.MCM'
1760 include 'COMMON.MCE'
1761 include 'COMMON.IOUNITS'
1763 character*320 mcmcard
1764 call card_concat(mcmcard)
1765 call readi(mcmcard,'MAXACC',maxacc,100)
1766 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1767 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1768 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1769 call readi(mcmcard,'MAXREPM',maxrepm,200)
1770 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1771 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1772 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1773 call reada(mcmcard,'E_UP',e_up,5.0D0)
1774 call reada(mcmcard,'DELTE',delte,0.1D0)
1775 call readi(mcmcard,'NSWEEP',nsweep,5)
1776 call readi(mcmcard,'NSTEPH',nsteph,0)
1777 call readi(mcmcard,'NSTEPC',nstepc,0)
1778 call reada(mcmcard,'TMIN',tmin,298.0D0)
1779 call reada(mcmcard,'TMAX',tmax,298.0D0)
1780 call readi(mcmcard,'NWINDOW',nwindow,0)
1781 call readi(mcmcard,'PRINT_MC',print_mc,0)
1782 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1783 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1784 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1785 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1786 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1787 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1788 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1789 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1790 if (nwindow.gt.0) then
1791 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1793 winlen(i)=winend(i)-winstart(i)+1
1796 if (tmax.lt.tmin) tmax=tmin
1797 if (tmax.eq.tmin) then
1801 if (nstepc.gt.0 .and. nsteph.gt.0) then
1802 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1803 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1805 C Probabilities of different move types
1806 sumpro_type(0)=0.0D0
1807 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1808 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1809 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1810 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1811 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1812 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1813 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1815 print *,'i',i,' sumprotype',sumpro_type(i)
1816 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1817 print *,'i',i,' sumprotype',sumpro_type(i)
1821 c----------------------------------------------------------------------------
1822 subroutine read_minim
1823 implicit real*8 (a-h,o-z)
1824 include 'DIMENSIONS'
1825 include 'COMMON.MINIM'
1826 include 'COMMON.IOUNITS'
1828 character*320 minimcard
1829 call card_concat(minimcard)
1830 call readi(minimcard,'MAXMIN',maxmin,2000)
1831 call readi(minimcard,'MAXFUN',maxfun,5000)
1832 call readi(minimcard,'MINMIN',minmin,maxmin)
1833 call readi(minimcard,'MINFUN',minfun,maxmin)
1834 call reada(minimcard,'TOLF',tolf,1.0D-2)
1835 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1836 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1837 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1838 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1839 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1840 & 'Options in energy minimization:'
1841 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1842 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1843 & 'MinMin:',MinMin,' MinFun:',MinFun,
1844 & ' TolF:',TolF,' RTolF:',RTolF
1847 c----------------------------------------------------------------------------
1848 subroutine read_angles(kanal,*)
1849 implicit real*8 (a-h,o-z)
1850 include 'DIMENSIONS'
1851 include 'COMMON.GEO'
1852 include 'COMMON.VAR'
1853 include 'COMMON.CHAIN'
1854 include 'COMMON.IOUNITS'
1855 include 'COMMON.CONTROL'
1856 c Read angles from input
1858 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1859 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1860 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1861 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1864 c 9/7/01 avoid 180 deg valence angle
1865 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1867 theta(i)=deg2rad*theta(i)
1868 phi(i)=deg2rad*phi(i)
1869 alph(i)=deg2rad*alph(i)
1870 omeg(i)=deg2rad*omeg(i)
1875 c----------------------------------------------------------------------------
1876 subroutine reada(rekord,lancuch,wartosc,default)
1878 character*(*) rekord,lancuch
1879 double precision wartosc,default
1882 iread=index(rekord,lancuch)
1883 if (iread.eq.0) then
1887 iread=iread+ilen(lancuch)+1
1888 read (rekord(iread:),*,err=10,end=10) wartosc
1893 c----------------------------------------------------------------------------
1894 subroutine readi(rekord,lancuch,wartosc,default)
1896 character*(*) rekord,lancuch
1897 integer wartosc,default
1900 iread=index(rekord,lancuch)
1901 if (iread.eq.0) then
1905 iread=iread+ilen(lancuch)+1
1906 read (rekord(iread:),*,err=10,end=10) wartosc
1911 c----------------------------------------------------------------------------
1912 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1915 integer tablica(dim),default
1916 character*(*) rekord,lancuch
1923 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1924 if (iread.eq.0) return
1925 iread=iread+ilen(lancuch)+1
1926 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1929 c----------------------------------------------------------------------------
1930 subroutine multreada(rekord,lancuch,tablica,dim,default)
1933 double precision tablica(dim),default
1934 character*(*) rekord,lancuch
1941 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1942 if (iread.eq.0) return
1943 iread=iread+ilen(lancuch)+1
1944 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1947 c----------------------------------------------------------------------------
1948 subroutine openunits
1949 implicit real*8 (a-h,o-z)
1950 include 'DIMENSIONS'
1953 character*16 form,nodename
1956 include 'COMMON.SETUP'
1957 include 'COMMON.IOUNITS'
1959 include 'COMMON.CONTROL'
1960 integer lenpre,lenpot,ilen,lentmp
1962 character*3 out1file_text,ucase
1967 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1968 call getenv_loc("PREFIX",prefix)
1970 call getenv_loc("POT",pot)
1971 call getenv_loc("DIRTMP",tmpdir)
1972 call getenv_loc("CURDIR",curdir)
1973 call getenv_loc("OUT1FILE",out1file_text)
1974 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1975 out1file_text=ucase(out1file_text)
1976 if (out1file_text(1:1).eq."Y") then
1979 out1file=fg_rank.gt.0
1984 if (lentmp.gt.0) then
1985 write (*,'(80(1h!))')
1986 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1987 write (*,'(80(1h!))')
1988 write (*,*)"All output files will be on node /tmp directory."
1990 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1991 if (me.eq.king) then
1992 write (*,*) "The master node is ",nodename
1993 else if (fg_rank.eq.0) then
1994 write (*,*) "I am the CG slave node ",nodename
1996 write (*,*) "I am the FG slave node ",nodename
1999 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2000 lenpre = lentmp+lenpre+1
2002 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2003 C Get the names and open the input files
2004 #if defined(WINIFL) || defined(WINPGI)
2005 open(1,file=pref_orig(:ilen(pref_orig))//
2006 & '.inp',status='old',readonly,shared)
2007 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2008 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2009 C Get parameter filenames and open the parameter files.
2010 call getenv_loc('BONDPAR',bondname)
2011 open (ibond,file=bondname,status='old',readonly,shared)
2012 call getenv_loc('THETPAR',thetname)
2013 open (ithep,file=thetname,status='old',readonly,shared)
2014 call getenv_loc('ROTPAR',rotname)
2015 open (irotam,file=rotname,status='old',readonly,shared)
2016 call getenv_loc('TORPAR',torname)
2017 open (itorp,file=torname,status='old',readonly,shared)
2018 call getenv_loc('TORDPAR',tordname)
2019 open (itordp,file=tordname,status='old',readonly,shared)
2020 call getenv_loc('FOURIER',fouriername)
2021 open (ifourier,file=fouriername,status='old',readonly,shared)
2022 call getenv_loc('ELEPAR',elename)
2023 open (ielep,file=elename,status='old',readonly,shared)
2024 call getenv_loc('SIDEPAR',sidename)
2025 open (isidep,file=sidename,status='old',readonly,shared)
2026 call getenv_loc('LIPTRANPAR',liptranname)
2027 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2028 #elif (defined CRAY) || (defined AIX)
2029 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2031 c print *,"Processor",myrank," opened file 1"
2032 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2033 c print *,"Processor",myrank," opened file 9"
2034 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2035 C Get parameter filenames and open the parameter files.
2036 call getenv_loc('BONDPAR',bondname)
2037 open (ibond,file=bondname,status='old',action='read')
2038 c print *,"Processor",myrank," opened file IBOND"
2039 call getenv_loc('THETPAR',thetname)
2040 open (ithep,file=thetname,status='old',action='read')
2042 call getenv_loc('THETPARPDB',thetname_pdb)
2043 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2045 c print *,"Processor",myrank," opened file ITHEP"
2046 call getenv_loc('ROTPAR',rotname)
2047 open (irotam,file=rotname,status='old',action='read')
2049 call getenv_loc('ROTPARPDB',rotname_pdb)
2050 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2052 c print *,"Processor",myrank," opened file IROTAM"
2053 call getenv_loc('TORPAR',torname)
2054 open (itorp,file=torname,status='old',action='read')
2055 c print *,"Processor",myrank," opened file ITORP"
2056 call getenv_loc('TORDPAR',tordname)
2057 open (itordp,file=tordname,status='old',action='read')
2058 c print *,"Processor",myrank," opened file ITORDP"
2059 call getenv_loc('SCCORPAR',sccorname)
2060 open (isccor,file=sccorname,status='old',action='read')
2061 c print *,"Processor",myrank," opened file ISCCOR"
2062 call getenv_loc('FOURIER',fouriername)
2063 open (ifourier,file=fouriername,status='old',action='read')
2064 c print *,"Processor",myrank," opened file IFOURIER"
2065 call getenv_loc('ELEPAR',elename)
2066 open (ielep,file=elename,status='old',action='read')
2067 c print *,"Processor",myrank," opened file IELEP"
2068 call getenv_loc('SIDEPAR',sidename)
2069 open (isidep,file=sidename,status='old',action='read')
2070 call getenv_loc('LIPTRANPAR',liptranname)
2071 open (iliptranpar,file=liptranname,status='old',action='read')
2072 c print *,"Processor",myrank," opened file ISIDEP"
2073 c print *,"Processor",myrank," opened parameter files"
2075 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2076 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2077 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2078 C Get parameter filenames and open the parameter files.
2079 call getenv_loc('BONDPAR',bondname)
2080 open (ibond,file=bondname,status='old')
2081 call getenv_loc('THETPAR',thetname)
2082 open (ithep,file=thetname,status='old')
2084 call getenv_loc('THETPARPDB',thetname_pdb)
2085 open (ithep_pdb,file=thetname_pdb,status='old')
2087 call getenv_loc('ROTPAR',rotname)
2088 open (irotam,file=rotname,status='old')
2090 call getenv_loc('ROTPARPDB',rotname_pdb)
2091 open (irotam_pdb,file=rotname_pdb,status='old')
2093 call getenv_loc('TORPAR',torname)
2094 open (itorp,file=torname,status='old')
2095 call getenv_loc('TORDPAR',tordname)
2096 open (itordp,file=tordname,status='old')
2097 call getenv_loc('SCCORPAR',sccorname)
2098 open (isccor,file=sccorname,status='old')
2099 call getenv_loc('FOURIER',fouriername)
2100 open (ifourier,file=fouriername,status='old')
2101 call getenv_loc('ELEPAR',elename)
2102 open (ielep,file=elename,status='old')
2103 call getenv_loc('SIDEPAR',sidename)
2104 open (isidep,file=sidename,status='old')
2105 call getenv_loc('LIPTRANPAR',liptranname)
2106 open (iliptranpar,file=liptranname,status='old')
2108 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2110 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2111 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2112 C Get parameter filenames and open the parameter files.
2113 call getenv_loc('BONDPAR',bondname)
2114 open (ibond,file=bondname,status='old',readonly)
2115 call getenv_loc('THETPAR',thetname)
2116 open (ithep,file=thetname,status='old',readonly)
2117 call getenv_loc('ROTPAR',rotname)
2118 open (irotam,file=rotname,status='old',readonly)
2119 call getenv_loc('TORPAR',torname)
2120 open (itorp,file=torname,status='old',readonly)
2121 call getenv_loc('TORDPAR',tordname)
2122 open (itordp,file=tordname,status='old',readonly)
2123 call getenv_loc('SCCORPAR',sccorname)
2124 open (isccor,file=sccorname,status='old',readonly)
2126 call getenv_loc('THETPARPDB',thetname_pdb)
2127 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2129 call getenv_loc('FOURIER',fouriername)
2130 open (ifourier,file=fouriername,status='old',readonly)
2131 call getenv_loc('ELEPAR',elename)
2132 open (ielep,file=elename,status='old',readonly)
2133 call getenv_loc('SIDEPAR',sidename)
2134 open (isidep,file=sidename,status='old',readonly)
2135 call getenv_loc('LIPTRANPAR',liptranname)
2136 open (iliptranpar,file=liptranname,status='old',action='read')
2138 call getenv_loc('ROTPARPDB',rotname_pdb)
2139 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2143 call getenv_loc('TUBEPAR',tubename)
2144 #if defined(WINIFL) || defined(WINPGI)
2145 open (itube,file=tubename,status='old',readonly,shared)
2146 #elif (defined CRAY) || (defined AIX)
2147 open (itube,file=tubename,status='old',action='read')
2149 open (itube,file=tubename,status='old')
2151 open (itube,file=tubename,status='old',readonly)
2156 C 8/9/01 In the newest version SCp interaction constants are read from a file
2157 C Use -DOLDSCP to use hard-coded constants instead.
2159 call getenv_loc('SCPPAR',scpname)
2160 #if defined(WINIFL) || defined(WINPGI)
2161 open (iscpp,file=scpname,status='old',readonly,shared)
2162 #elif (defined CRAY) || (defined AIX)
2163 open (iscpp,file=scpname,status='old',action='read')
2165 open (iscpp,file=scpname,status='old')
2167 open (iscpp,file=scpname,status='old',readonly)
2170 call getenv_loc('PATTERN',patname)
2171 #if defined(WINIFL) || defined(WINPGI)
2172 open (icbase,file=patname,status='old',readonly,shared)
2173 #elif (defined CRAY) || (defined AIX)
2174 open (icbase,file=patname,status='old',action='read')
2176 open (icbase,file=patname,status='old')
2178 open (icbase,file=patname,status='old',readonly)
2181 C Open output file only for CG processes
2182 c print *,"Processor",myrank," fg_rank",fg_rank
2183 if (fg_rank.eq.0) then
2185 if (nodes.eq.1) then
2188 npos = dlog10(dfloat(nodes-1))+1
2190 if (npos.lt.3) npos=3
2191 write (liczba,'(i1)') npos
2192 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2194 write (liczba,form) me
2195 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2196 & liczba(:ilen(liczba))
2197 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2199 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2201 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2202 & liczba(:ilen(liczba))//'.mol2'
2203 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2204 & liczba(:ilen(liczba))//'.stat'
2206 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2207 & //liczba(:ilen(liczba))//'.stat')
2208 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2211 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2212 & liczba(:ilen(liczba))//'.const'
2217 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2218 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2219 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2220 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2221 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2223 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2225 rest2name=prefix(:ilen(prefix))//'.rst'
2227 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2230 #if defined(AIX) || defined(PGI) || defined(CRAY)
2231 if (me.eq.king .or. .not. out1file) then
2232 open(iout,file=outname,status='unknown')
2234 open(iout,file="/dev/null",status="unknown")
2238 if (fg_rank.gt.0) then
2239 write (liczba,'(i3.3)') myrank/nfgtasks
2240 write (ll,'(bz,i3.3)') fg_rank
2241 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2247 open(igeom,file=intname,status='unknown',position='append')
2248 open(ipdb,file=pdbname,status='unknown')
2249 open(imol2,file=mol2name,status='unknown')
2250 open(istat,file=statname,status='unknown',position='append')
2252 c1out open(iout,file=outname,status='unknown')
2255 if (me.eq.king .or. .not.out1file)
2256 & open(iout,file=outname,status='unknown')
2258 if (fg_rank.gt.0) then
2259 write (liczba,'(i3.3)') myrank/nfgtasks
2260 write (ll,'(bz,i3.3)') fg_rank
2261 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2266 open(igeom,file=intname,status='unknown',access='append')
2267 open(ipdb,file=pdbname,status='unknown')
2268 open(imol2,file=mol2name,status='unknown')
2269 open(istat,file=statname,status='unknown',access='append')
2271 c1out open(iout,file=outname,status='unknown')
2274 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2275 csa_seed=prefix(:lenpre)//'.CSA.seed'
2276 csa_history=prefix(:lenpre)//'.CSA.history'
2277 csa_bank=prefix(:lenpre)//'.CSA.bank'
2278 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2279 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2280 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2281 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2282 csa_int=prefix(:lenpre)//'.int'
2283 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2284 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2285 csa_in=prefix(:lenpre)//'.CSA.in'
2286 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2289 write (iout,'(80(1h-))')
2290 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2291 write (iout,'(80(1h-))')
2292 write (iout,*) "Input file : ",
2293 & pref_orig(:ilen(pref_orig))//'.inp'
2294 write (iout,*) "Output file : ",
2295 & outname(:ilen(outname))
2297 write (iout,*) "Sidechain potential file : ",
2298 & sidename(:ilen(sidename))
2300 write (iout,*) "SCp potential file : ",
2301 & scpname(:ilen(scpname))
2303 write (iout,*) "Electrostatic potential file : ",
2304 & elename(:ilen(elename))
2305 write (iout,*) "Cumulant coefficient file : ",
2306 & fouriername(:ilen(fouriername))
2307 write (iout,*) "Torsional parameter file : ",
2308 & torname(:ilen(torname))
2309 write (iout,*) "Double torsional parameter file : ",
2310 & tordname(:ilen(tordname))
2311 write (iout,*) "SCCOR parameter file : ",
2312 & sccorname(:ilen(sccorname))
2313 write (iout,*) "Bond & inertia constant file : ",
2314 & bondname(:ilen(bondname))
2315 write (iout,*) "Bending-torsion parameter file : ",
2316 & thetname(:ilen(thetname))
2317 write (iout,*) "Rotamer parameter file : ",
2318 & rotname(:ilen(rotname))
2319 write (iout,*) "Thetpdb parameter file : ",
2320 & thetname_pdb(:ilen(thetname_pdb))
2321 write (iout,*) "Threading database : ",
2322 & patname(:ilen(patname))
2324 &write (iout,*)" DIRTMP : ",
2326 write (iout,'(80(1h-))')
2330 c----------------------------------------------------------------------------
2331 subroutine card_concat(card)
2332 implicit real*8 (a-h,o-z)
2333 include 'DIMENSIONS'
2334 include 'COMMON.IOUNITS'
2336 character*80 karta,ucase
2338 read (inp,'(a)') karta
2341 do while (karta(80:80).eq.'&')
2342 card=card(:ilen(card)+1)//karta(:79)
2343 read (inp,'(a)') karta
2346 card=card(:ilen(card)+1)//karta
2349 c----------------------------------------------------------------------------------
2351 implicit real*8 (a-h,o-z)
2352 include 'DIMENSIONS'
2353 include 'COMMON.CHAIN'
2354 include 'COMMON.IOUNITS'
2356 open(irest2,file=rest2name,status='unknown')
2357 read(irest2,*) totT,EK,potE,totE,t_bath
2360 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2363 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2366 read (irest2,*) iset
2371 c---------------------------------------------------------------------------------
2372 subroutine read_fragments
2373 implicit real*8 (a-h,o-z)
2374 include 'DIMENSIONS'
2378 include 'COMMON.SETUP'
2379 include 'COMMON.CHAIN'
2380 include 'COMMON.IOUNITS'
2382 include 'COMMON.CONTROL'
2383 read(inp,*) nset,nfrag,npair,nfrag_back
2384 loc_qlike=(nfrag_back.lt.0)
2385 nfrag_back=iabs(nfrag_back)
2386 if(me.eq.king.or..not.out1file)
2387 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2388 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2390 read(inp,*) mset(iset)
2392 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2394 if(me.eq.king.or..not.out1file)
2395 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2396 & ifrag(2,i,iset), qinfrag(i,iset)
2399 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2401 if(me.eq.king.or..not.out1file)
2402 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2403 & ipair(2,i,iset), qinpair(i,iset)
2407 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2408 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2409 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2410 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2411 if(me.eq.king.or..not.out1file)
2412 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2413 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2414 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2415 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2419 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2420 & wfrag_back(3,i,iset),
2421 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2422 if(me.eq.king.or..not.out1file)
2423 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2424 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2430 C---------------------------------------------------------------------------
2431 subroutine read_afminp
2432 implicit real*8 (a-h,o-z)
2433 include 'DIMENSIONS'
2437 include 'COMMON.SETUP'
2438 include 'COMMON.CONTROL'
2439 include 'COMMON.CHAIN'
2440 include 'COMMON.IOUNITS'
2441 include 'COMMON.SBRIDGE'
2442 character*320 afmcard
2444 call card_concat(afmcard)
2445 call readi(afmcard,"BEG",afmbeg,0)
2446 call readi(afmcard,"END",afmend,0)
2447 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2448 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2449 print *,'FORCE=' ,forceAFMconst
2450 CCCC NOW PROPERTIES FOR AFM
2453 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2455 distafminit=dsqrt(distafminit)
2456 print *,'initdist',distafminit
2459 c-------------------------------------------------------------------------------
2460 subroutine read_saxs_constr
2461 implicit real*8 (a-h,o-z)
2462 include 'DIMENSIONS'
2466 include 'COMMON.SETUP'
2467 include 'COMMON.CONTROL'
2468 include 'COMMON.CHAIN'
2469 include 'COMMON.IOUNITS'
2470 include 'COMMON.SBRIDGE'
2471 double precision cm(3)
2473 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2475 if (saxs_mode.eq.0) then
2476 c SAXS distance distribution
2478 read(inp,*) distsaxs(i),Psaxs(i)
2482 Cnorm = Cnorm + Psaxs(i)
2484 write (iout,*) "Cnorm",Cnorm
2486 Psaxs(i)=Psaxs(i)/Cnorm
2488 write (iout,*) "Normalized distance distribution from SAXS"
2490 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2494 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2496 write (iout,*) "Wsaxs0",Wsaxs0
2500 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2507 cm(j)=cm(j)+Csaxs(j,i)
2515 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2518 write (iout,*) "SAXS sphere coordinates"
2520 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2526 c-------------------------------------------------------------------------------
2527 subroutine read_dist_constr
2528 implicit real*8 (a-h,o-z)
2529 include 'DIMENSIONS'
2533 include 'COMMON.SETUP'
2534 include 'COMMON.CONTROL'
2535 include 'COMMON.CHAIN'
2536 include 'COMMON.IOUNITS'
2537 include 'COMMON.SBRIDGE'
2538 include 'COMMON.INTERACT'
2539 integer ifrag_(2,100),ipair_(2,1000)
2540 double precision wfrag_(100),wpair_(1000)
2541 character*5000 controlcard
2542 logical normalize,next
2544 double precision scal_bfac
2545 double precision xlink(4,0:4) /
2547 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2548 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2549 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2550 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2551 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2552 c print *, "WCHODZE"
2553 write (iout,*) "Calling read_dist_constr"
2554 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2556 restr_on_coord=.false.
2565 call card_concat(controlcard)
2566 next = index(controlcard,"NEXT").gt.0
2567 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2568 write (iout,*) "restr_type",restr_type
2569 call readi(controlcard,"NFRAG",nfrag_,0)
2570 call readi(controlcard,"NFRAG",nfrag_,0)
2571 call readi(controlcard,"NPAIR",npair_,0)
2572 call readi(controlcard,"NDIST",ndist_,0)
2573 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2574 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2575 if (restr_type.eq.10)
2576 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2577 if (restr_type.eq.12)
2578 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2579 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2580 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2581 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2582 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2583 normalize = index(controlcard,"NORMALIZE").gt.0
2584 write (iout,*) "WBOLTZD",wboltzd
2585 write (iout,*) "SCAL_PEAK",scal_peak
2586 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2587 write (iout,*) "IFRAG"
2589 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2591 write (iout,*) "IPAIR"
2593 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2595 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2597 & "Distance restraints as generated from reference structure"
2599 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2600 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2601 & ifrag_(2,i)=nstart_sup+nsup-1
2602 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2604 if (wfrag_(i).eq.0.0d0) cycle
2605 do j=ifrag_(1,i),ifrag_(2,i)-1
2606 do k=j+1,ifrag_(2,i)
2607 c write (iout,*) "j",j," k",k
2609 if (restr_type.eq.1) then
2615 forcon(nhpb)=wfrag_(i)
2616 else if (constr_dist.eq.2) then
2617 if (ddjk.le.dist_cut) then
2623 forcon(nhpb)=wfrag_(i)
2625 else if (restr_type.eq.3) then
2631 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2634 if (.not.out1file .or. me.eq.king)
2635 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2636 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2638 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2639 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2645 if (wpair_(i).eq.0.0d0) cycle
2653 do j=ifrag_(1,ii),ifrag_(2,ii)
2654 do k=ifrag_(1,jj),ifrag_(2,jj)
2656 if (restr_type.eq.1) then
2662 forcon(nhpb)=wpair_(i)
2663 else if (constr_dist.eq.2) then
2664 if (ddjk.le.dist_cut) then
2670 forcon(nhpb)=wpair_(i)
2672 else if (restr_type.eq.3) then
2678 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2681 if (.not.out1file .or. me.eq.king)
2682 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2683 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2685 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2686 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2693 write (iout,*) "Distance restraints as read from input"
2695 if (restr_type.eq.12) then
2696 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2697 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2698 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2699 & fordepth_peak(nhpb_peak+1),npeak
2700 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2701 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2702 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2703 c & fordepth_peak(nhpb_peak+1),npeak
2704 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2705 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2706 nhpb_peak=nhpb_peak+1
2707 irestr_type_peak(nhpb_peak)=12
2708 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2711 if (.not.out1file .or. me.eq.king)
2712 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2713 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2714 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2715 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2716 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2718 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2719 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2720 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2721 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2722 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2724 if (ibecarb_peak(nhpb_peak).eq.3) then
2725 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2726 else if (ibecarb_peak(nhpb_peak).eq.2) then
2727 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2728 else if (ibecarb_peak(nhpb_peak).eq.1) then
2729 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2730 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2732 else if (restr_type.eq.11) then
2733 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2734 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2735 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2736 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2738 irestr_type(nhpb)=11
2740 if (.not.out1file .or. me.eq.king)
2741 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2742 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2743 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2745 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2746 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2747 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2749 c if (ibecarb(nhpb).gt.0) then
2750 c ihpb(nhpb)=ihpb(nhpb)+nres
2751 c jhpb(nhpb)=jhpb(nhpb)+nres
2753 if (ibecarb(nhpb).eq.3) then
2754 ihpb(nhpb)=ihpb(nhpb)+nres
2755 else if (ibecarb(nhpb).eq.2) then
2756 ihpb(nhpb)=ihpb(nhpb)+nres
2757 else if (ibecarb(nhpb).eq.1) then
2758 ihpb(nhpb)=ihpb(nhpb)+nres
2759 jhpb(nhpb)=jhpb(nhpb)+nres
2761 else if (restr_type.eq.10) then
2762 c Cross-lonk Markov-like potential
2763 call card_concat(controlcard)
2764 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2765 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2767 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2768 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2769 if (index(controlcard,"ZL").gt.0) then
2771 else if (index(controlcard,"ADH").gt.0) then
2773 else if (index(controlcard,"PDH").gt.0) then
2775 else if (index(controlcard,"DSS").gt.0) then
2780 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2781 & xlink(1,link_type))
2782 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2783 & xlink(2,link_type))
2784 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2785 & xlink(3,link_type))
2786 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2787 & xlink(4,link_type))
2788 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2789 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2790 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2791 if (forcon(nhpb+1).le.0.0d0 .or.
2792 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2794 irestr_type(nhpb)=10
2795 if (ibecarb(nhpb).eq.3) then
2796 jhpb(nhpb)=jhpb(nhpb)+nres
2797 else if (ibecarb(nhpb).eq.2) then
2798 ihpb(nhpb)=ihpb(nhpb)+nres
2799 else if (ibecarb(nhpb).eq.1) then
2800 ihpb(nhpb)=ihpb(nhpb)+nres
2801 jhpb(nhpb)=jhpb(nhpb)+nres
2804 if (.not.out1file .or. me.eq.king)
2805 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2806 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2807 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2810 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2811 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2812 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2817 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2818 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2819 if (forcon(nhpb+1).gt.0.0d0) then
2821 if (dhpb1(nhpb).eq.0.0d0) then
2826 if (ibecarb(nhpb).eq.3) then
2827 jhpb(nhpb)=jhpb(nhpb)+nres
2828 else if (ibecarb(nhpb).eq.2) then
2829 ihpb(nhpb)=ihpb(nhpb)+nres
2830 else if (ibecarb(nhpb).eq.1) then
2831 ihpb(nhpb)=ihpb(nhpb)+nres
2832 jhpb(nhpb)=jhpb(nhpb)+nres
2834 if (dhpb(nhpb).eq.0.0d0)
2835 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2838 if (.not.out1file .or. me.eq.king)
2839 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2840 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2842 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2843 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2846 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2847 C if (forcon(nhpb+1).gt.0.0d0) then
2849 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2852 if (restr_type.eq.4) then
2853 write (iout,*) "The BFAC array"
2855 write (iout,'(i5,f10.5)') i,bfac(i)
2858 if (itype(i).eq.ntyp1) cycle
2860 if (itype(j).eq.ntyp1) cycle
2861 if (itype(i).eq.10) then
2866 if (itype(j).eq.10) then
2876 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2879 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2880 ihpb(nhpb)=i+nres*ii
2881 jhpb(nhpb)=j+nres*jj
2882 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2884 if (.not.out1file .or. me.eq.king) then
2885 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2886 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2887 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2891 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2892 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2893 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2903 if (restr_type.eq.5) then
2904 restr_on_coord=.true.
2906 if (itype(i).eq.ntyp1) cycle
2907 bfac(i)=(scal_bfac/bfac(i))**2
2916 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2917 & fordepthmax=fordepth(i)
2920 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2925 c-------------------------------------------------------------------------------
2927 subroutine flush(iu)
2932 subroutine flush(iu)
2937 c------------------------------------------------------------------------------
2938 subroutine copy_to_tmp(source)
2939 include "DIMENSIONS"
2940 include "COMMON.IOUNITS"
2941 character*(*) source
2942 character* 256 tmpfile
2946 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2947 inquire(file=tmpfile,exist=ex)
2949 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2950 & " to temporary directory..."
2951 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2952 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2956 c------------------------------------------------------------------------------
2957 subroutine move_from_tmp(source)
2958 include "DIMENSIONS"
2959 include "COMMON.IOUNITS"
2960 character*(*) source
2963 write (*,*) "Moving ",source(:ilen(source)),
2964 & " from temporary directory to working directory"
2965 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2966 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2969 c------------------------------------------------------------------------------
2970 subroutine random_init(seed)
2972 C Initialize random number generator
2974 implicit real*8 (a-h,o-z)
2975 include 'DIMENSIONS'
2978 logical OKRandom, prng_restart
2980 integer iseed_array(4)
2982 include 'COMMON.IOUNITS'
2983 include 'COMMON.TIME1'
2984 include 'COMMON.THREAD'
2985 include 'COMMON.SBRIDGE'
2986 include 'COMMON.CONTROL'
2987 include 'COMMON.MCM'
2988 include 'COMMON.MAP'
2989 include 'COMMON.HEADER'
2990 include 'COMMON.CSA'
2991 include 'COMMON.CHAIN'
2992 include 'COMMON.MUCA'
2994 include 'COMMON.FFIELD'
2995 include 'COMMON.SETUP'
2996 iseed=-dint(dabs(seed))
2997 if (iseed.eq.0) then
2998 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2999 & 'Random seed undefined. The program will stop.'
3000 write (*,'(/80(1h*)/20x,a/80(1h*))')
3001 & 'Random seed undefined. The program will stop.'
3003 call mpi_finalize(mpi_comm_world,ierr)
3005 stop 'Bad random seed.'
3008 if (fg_rank.eq.0) then
3012 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3013 OKRandom = prng_restart(me,iseed)
3016 tmp=65536.0d0**(4-i)
3017 iseed_array(i) = dint(seed/tmp)
3018 seed=seed-iseed_array(i)*tmp
3021 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3022 & (iseed_array(i),i=1,4)
3023 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3024 & (iseed_array(i),i=1,4)
3025 OKRandom = prng_restart(me,iseed_array)
3028 c r1 = prng_next(me)
3029 r1=ran_number(0.0D0,1.0D0)
3031 & write (iout,*) 'ran_num',r1
3032 if (r1.lt.0.0d0) OKRandom=.false.
3034 if (.not.OKRandom) then
3035 write (iout,*) 'PRNG IS NOT WORKING!!!'
3036 print *,'PRNG IS NOT WORKING!!!'
3039 call mpi_abort(mpi_comm_world,error_msg,ierr)
3042 write (iout,*) 'too many processors for parallel prng'
3043 write (*,*) 'too many processors for parallel prng'
3051 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)